Index: branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90
===================================================================
--- branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90	(revision 8873)
+++ branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90	(revision 8882)
@@ -5,10 +5,10 @@
    !!                 turbulent closure parameterization
    !!======================================================================
-   !! History :   3.0  !  2009-09  (G. Reffray)  Original code
-   !!             3.3  !  2010-10  (C. Bricaud)  Add in the reference
+   !! History :  3.0  !  2009-09  (G. Reffray)  Original code
+   !!            3.3  !  2010-10  (C. Bricaud)  Add in the reference
+   !!            4.0  !  2017-04  (G. Madec)  remove CPP keys & avm at t-point only 
+   !!             -   !  2017-05  (G. Madec)  add top friction as boundary condition
    !!----------------------------------------------------------------------
-#if defined key_zdfgls
-   !!----------------------------------------------------------------------
-   !!   'key_zdfgls'                 Generic Length Scale vertical physics
+
    !!----------------------------------------------------------------------
    !!   zdf_gls       : update momentum and tracer Kz from a gls scheme
@@ -20,13 +20,13 @@
    USE domvvl         ! ocean space and time domain : variable volume layer
    USE zdf_oce        ! ocean vertical physics
-   USE zdfbfr         ! bottom friction (only for rn_bfrz0)
+   USE zdfdrg  , ONLY : r_z0_top , r_z0_bot   ! top/bottom roughness
+   USE zdfdrg  , ONLY : rCdU_top , rCdU_bot   ! top/bottom friction
    USE sbc_oce        ! surface boundary condition: ocean
    USE phycst         ! physical constants
    USE zdfmxl         ! mixed layer
-   USE sbcwave ,  ONLY: hsw   ! significant wave height
+   USE sbcwave , ONLY : hsw   ! significant wave height
    !
    USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
    USE lib_mpp        ! MPP manager
-   USE wrk_nemo       ! work arrays
    USE prtctl         ! Print control
    USE in_out_manager ! I/O manager
@@ -38,14 +38,14 @@
    PRIVATE
 
-   PUBLIC   zdf_gls        ! routine called in step module
-   PUBLIC   zdf_gls_init   ! routine called in opa module
-   PUBLIC   gls_rst        ! routine called in step module
-
-   LOGICAL , PUBLIC, PARAMETER ::   lk_zdfgls = .TRUE.   !: TKE vertical mixing flag
+   PUBLIC   zdf_gls        ! called in zdfphy
+   PUBLIC   zdf_gls_init   ! called in zdfphy
+   PUBLIC   gls_rst        ! called in zdfphy
+
    !
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   mxln    !: now mixing length
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hmxl_n    !: now mixing length
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   zwall   !: wall function
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ustars2 !: Squared surface velocity scale at T-points
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ustarb2 !: Squared bottom  velocity scale at T-points
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ustar2_surf !: Squared surface velocity scale at T-points
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ustar2_top  !: Squared top     velocity scale at T-points
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ustar2_bot  !: Squared bottom  velocity scale at T-points
 
    !                              !! ** Namelist  namzdf_gls  **
@@ -102,4 +102,6 @@
    REAL(wp) ::   rsc_tke, rsc_psi, rpsi1, rpsi2, rpsi3, rsc_psi0  !     -           -           -        -
    REAL(wp) ::   rpsi3m, rpsi3p, rpp, rmm, rnn                    !     -           -           -        -
+   !
+   REAL(wp) ::   r2_3 = 2._wp/3._wp   ! constant=2/3
 
    !! * Substitutions
@@ -116,6 +118,6 @@
       !!                ***  FUNCTION zdf_gls_alloc  ***
       !!----------------------------------------------------------------------
-      ALLOCATE( mxln(jpi,jpj,jpk), zwall(jpi,jpj,jpk) ,     &
-         &      ustars2(jpi,jpj) , ustarb2(jpi,jpj)   , STAT= zdf_gls_alloc )
+      ALLOCATE( hmxl_n(jpi,jpj,jpk) , ustar2_surf(jpi,jpj) ,                     &
+         &      zwall (jpi,jpj,jpk) , ustar2_top (jpi,jpj) , ustar2_bot(jpi,jpj) , STAT= zdf_gls_alloc )
          !
       IF( lk_mpp             )   CALL mpp_sum ( zdf_gls_alloc )
@@ -124,5 +126,5 @@
 
 
-   SUBROUTINE zdf_gls( kt )
+   SUBROUTINE zdf_gls( kt, p_sh2, p_avm, p_avt )
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE zdf_gls  ***
@@ -131,94 +133,90 @@
       !!              coefficients using the GLS turbulent closure scheme.
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt ! ocean time step
-      INTEGER  ::   ji, jj, jk, ibot, ibotm1, dir  ! dummy loop arguments
-      REAL(wp) ::   zesh2, zsigpsi, zcoef, zex1, zex2   ! local scalars
-      REAL(wp) ::   ztx2, zty2, zup, zdown, zcof        !   -      - 
-      REAL(wp) ::   zratio, zrn2, zflxb, sh             !   -      -
+      INTEGER                   , INTENT(in   ) ::   kt             ! ocean time step
+      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   p_sh2          ! shear production term
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   p_avm, p_avt   !  momentum and tracer Kz (w-points)
+      !
+      INTEGER  ::   ji, jj, jk    ! dummy loop arguments
+      INTEGER  ::   ibot, ibotm1  ! local integers
+      INTEGER  ::   itop, itopp1  !   -       -
+      REAL(wp) ::   zesh2, zsigpsi, zcoef, zex1 , zex2  ! local scalars
+      REAL(wp) ::   ztx2, zty2, zup, zdown, zcof, zdir  !   -      - 
+      REAL(wp) ::   zratio, zrn2, zflxb, sh     , z_en  !   -      -
       REAL(wp) ::   prod, buoy, diss, zdiss, sm         !   -      -
-      REAL(wp) ::   gh, gm, shr, dif, zsqen, zav        !   -      -
-      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zdep
-      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zkar
-      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zflxs       ! Turbulence fluxed induced by internal waves 
-      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zhsro       ! Surface roughness (surface waves)
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   eb          ! tke at time before
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   mxlb        ! mixing length at time before
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   shear       ! vertical shear
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   eps         ! dissipation rate
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zwall_psi   ! Wall function use in the wb case (ln_sigpsi)
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   psi         ! psi at time now
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   z_elem_a    ! element of the first  matrix diagonal
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   z_elem_b    ! element of the second matrix diagonal
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   z_elem_c    ! element of the third  matrix diagonal
+      REAL(wp) ::   gh, gm, shr, dif, zsqen, zavt, zavm !   -      -
+      REAL(wp) ::   zmsku, zmskv                        !   -      -
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zdep
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zkar
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zflxs       ! Turbulence fluxed induced by internal waves 
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zhsro       ! Surface roughness (surface waves)
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   eb          ! tke at time before
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   hmxl_b      ! mixing length at time before
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   eps         ! dissipation rate
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwall_psi   ! Wall function use in the wb case (ln_sigpsi)
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   psi         ! psi at time now
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zd_lw, zd_up, zdiag   ! lower, upper  and diagonal of the matrix
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zstt, zstm  ! stability function on tracer and momentum
       !!--------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )  CALL timing_start('zdf_gls')
-      !
-      CALL wrk_alloc( jpi,jpj,       zdep, zkar, zflxs, zhsro )
-      CALL wrk_alloc( jpi,jpj,jpk,   eb, mxlb, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi  )
-      
+      IF( ln_timing )   CALL timing_start('zdf_gls')
+      !
       ! Preliminary computing
 
-      ustars2 = 0._wp   ;   ustarb2 = 0._wp   ;   psi  = 0._wp   ;   zwall_psi = 0._wp
-
-      IF( kt /= nit000 ) THEN   ! restore before value to compute tke
-         avt (:,:,:) = avt_k (:,:,:)
-         avm (:,:,:) = avm_k (:,:,:)
-         avmu(:,:,:) = avmu_k(:,:,:)
-         avmv(:,:,:) = avmv_k(:,:,:) 
-      ENDIF
-
-      ! Compute surface and bottom friction at T-points
+      ustar2_surf(:,:) = 0._wp   ;         psi(:,:,:) = 0._wp   
+      ustar2_top (:,:) = 0._wp   ;   zwall_psi(:,:,:) = 0._wp
+      ustar2_bot (:,:) = 0._wp
+
+      ! Compute surface, top and bottom friction at T-points
       DO jj = 2, jpjm1          
          DO ji = fs_2, fs_jpim1   ! vector opt.         
             !
             ! surface friction
-            ustars2(ji,jj) = r1_rau0 * taum(ji,jj) * tmask(ji,jj,1)
+            ustar2_surf(ji,jj) = r1_rau0 * taum(ji,jj) * tmask(ji,jj,1)
             !   
-            ! bottom friction (explicit before friction)        
-            ! Note that we chose here not to bound the friction as in dynbfr)   
-            ztx2 = (  bfrua(ji,jj)  * ub(ji,jj,mbku(ji,jj)) + bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj))  )   &         
-               & * ( 1._wp - 0.5_wp * umask(ji,jj,1) * umask(ji-1,jj,1)  )      
-            zty2 = (  bfrva(ji,jj)  * vb(ji,jj,mbkv(ji,jj)) + bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1))  )   &         
-               & * ( 1._wp - 0.5_wp * vmask(ji,jj,1) * vmask(ji,jj-1,1)  )      
-            ustarb2(ji,jj) = SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)         
-         END DO         
-      END DO    
-
-      ! Set surface roughness length
-      SELECT CASE ( nn_z0_met )
-      !
-      CASE ( 0 )             ! Constant roughness          
+!!gm Rq we may add here r_ke0(_top/_bot) ?  ==>> think about that...
+          ! bottom friction (explicit before friction)
+          zmsku = ( 2._wp - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) )
+          zmskv = ( 2._wp - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) )     ! (CAUTION: CdU<0)
+          ustar2_bot(ji,jj) = - rCdU_bot(ji,jj) * SQRT(  ( zmsku*( ub(ji,jj,mbkt(ji,jj))+ub(ji-1,jj,mbkt(ji,jj)) ) )**2  &
+             &                                         + ( zmskv*( vb(ji,jj,mbkt(ji,jj))+vb(ji,jj-1,mbkt(ji,jj)) ) )**2  )
+         END DO
+      END DO
+      IF( ln_isfcav ) THEN       !top friction
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) )
+               zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) )     ! (CAUTION: CdU<0)
+               ustar2_top(ji,jj) = - rCdU_top(ji,jj) * SQRT(  ( zmsku*( ub(ji,jj,mikt(ji,jj))+ub(ji-1,jj,mikt(ji,jj)) ) )**2  &
+                  &                                         + ( zmskv*( vb(ji,jj,mikt(ji,jj))+vb(ji,jj-1,mikt(ji,jj)) ) )**2  )
+            END DO
+         END DO
+      ENDIF
+   
+      SELECT CASE ( nn_z0_met )      !==  Set surface roughness length  ==!
+      CASE ( 0 )                          ! Constant roughness          
          zhsro(:,:) = rn_hsro
       CASE ( 1 )             ! Standard Charnock formula
-         zhsro(:,:) = MAX(rsbc_zs1 * ustars2(:,:), rn_hsro)
+         zhsro(:,:) = MAX( rsbc_zs1 * ustar2_surf(:,:) , rn_hsro )
       CASE ( 2 )             ! Roughness formulae according to Rascle et al., Ocean Modelling (2008)
-         zdep(:,:)  = 30.*TANH(2.*0.3/(28.*SQRT(MAX(ustars2(:,:),rsmall))))             ! Wave age (eq. 10)
-         zhsro(:,:) = MAX(rsbc_zs2 * ustars2(:,:) * zdep(:,:)**1.5, rn_hsro) ! zhsro = rn_frac_hs * Hsw (eq. 11)
+!!gm faster coding : the 2 comment lines should be used
+!!gm         zcof = 2._wp * 0.6_wp / 28._wp
+!!gm         zdep(:,:)  = 30._wp * TANH(  zcof/ SQRT( MAX(ustar2_surf(:,:),rsmall) )  )       ! Wave age (eq. 10)
+         zdep (:,:) = 30.*TANH( 2.*0.3/(28.*SQRT(MAX(ustar2_surf(:,:),rsmall))) )         ! Wave age (eq. 10)
+         zhsro(:,:) = MAX(rsbc_zs2 * ustar2_surf(:,:) * zdep(:,:)**1.5, rn_hsro)          ! zhsro = rn_frac_hs * Hsw (eq. 11)
       CASE ( 3 )             ! Roughness given by the wave model (coupled or read in file)
-         zhsro(:,:) = hsw(:,:)
+         zhsro(:,:) = rn_frac_hs * hsw(:,:)   ! (rn_frac_hs=1.6 see Eq. (5) of Rascle et al. 2008 )
       END SELECT
-
-      ! Compute shear and dissipation rate
-      DO jk = 2, jpkm1
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               avmu(ji,jj,jk) = avmu(ji,jj,jk) * ( un(ji,jj,jk-1) - un(ji,jj,jk) )   &
-                  &                            * ( ub(ji,jj,jk-1) - ub(ji,jj,jk) )   &
-                  &                            / (  e3uw_n(ji,jj,jk) * e3uw_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) )   &
-                  &                            / (  e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) )
-               eps(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT(en(ji,jj,jk)) / mxln(ji,jj,jk)
+      !
+      DO jk = 2, jpkm1              !==  Compute dissipation rate  ==!
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1
+               eps(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / hmxl_n(ji,jj,jk)
             END DO
          END DO
       END DO
-      !
-      ! Lateral boundary conditions (avmu,avmv) (sign unchanged)
-      CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )
 
       ! Save tke at before time step
-      eb  (:,:,:) = en  (:,:,:)
-      mxlb(:,:,:) = mxln(:,:,:)
+      eb    (:,:,:) = en    (:,:,:)
+      hmxl_b(:,:,:) = hmxl_n(:,:,:)
 
       IF( nn_clos == 0 ) THEN    ! Mellor-Yamada
@@ -226,5 +224,5 @@
             DO jj = 2, jpjm1 
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zup   = mxln(ji,jj,jk) * gdepw_n(ji,jj,mbkt(ji,jj)+1)
+                  zup   = hmxl_n(ji,jj,jk) * gdepw_n(ji,jj,mbkt(ji,jj)+1)
                   zdown = vkarmn * gdepw_n(ji,jj,jk) * ( -gdepw_n(ji,jj,jk) + gdepw_n(ji,jj,mbkt(ji,jj)+1) )
                   zcoef = ( zup / MAX( zdown, rsmall ) )
@@ -245,25 +243,23 @@
       ! The surface boundary condition are set after
       ! The bottom boundary condition are also set after. In standard e(bottom)=0.
-      ! z_elem_b : diagonal z_elem_c : upper diagonal z_elem_a : lower diagonal
+      ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal
       ! Warning : after this step, en : right hand side of the matrix
 
       DO jk = 2, jpkm1
          DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               !
-               ! shear prod. at w-point weightened by mask
-               shear(ji,jj,jk) =  ( avmu(ji-1,jj,jk) + avmu(ji,jj,jk) ) / MAX( 1.e0 , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   &
-                  &             + ( avmv(ji,jj-1,jk) + avmv(ji,jj,jk) ) / MAX( 1.e0 , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )
-               !
-               ! stratif. destruction
-               buoy = - avt(ji,jj,jk) * rn2(ji,jj,jk)
-               !
-               ! shear prod. - stratif. destruction
-               diss = eps(ji,jj,jk)
-               !
-               dir = 0.5_wp + SIGN( 0.5_wp, shear(ji,jj,jk) + buoy )   ! dir =1(=0) if shear(ji,jj,jk)+buoy >0(<0)
-               !
-               zesh2 = dir*(shear(ji,jj,jk)+buoy)+(1._wp-dir)*shear(ji,jj,jk)          ! production term
-               zdiss = dir*(diss/en(ji,jj,jk))   +(1._wp-dir)*(diss-buoy)/en(ji,jj,jk) ! dissipation term
+            DO ji = 2, jpim1
+               !
+               buoy = - p_avt(ji,jj,jk) * rn2(ji,jj,jk)     ! stratif. destruction
+               !
+               diss = eps(ji,jj,jk)                         ! dissipation
+               !
+               zdir = 0.5_wp + SIGN( 0.5_wp, p_sh2(ji,jj,jk) + buoy )   ! zdir =1(=0) if shear(ji,jj,jk)+buoy >0(<0)
+               !
+               zesh2 = zdir*(p_sh2(ji,jj,jk)+buoy)+(1._wp-zdir)*p_sh2(ji,jj,jk)          ! production term
+               zdiss = zdir*(diss/en(ji,jj,jk))   +(1._wp-zdir)*(diss-buoy)/en(ji,jj,jk) ! dissipation term
+!!gm better coding, identical results
+!               zesh2 =   p_sh2(ji,jj,jk) + zdir*buoy               ! production term
+!               zdiss = ( diss - (1._wp-zdir)*buoy ) / en(ji,jj,jk) ! dissipation term
+!!gm
                !
                ! Compute a wall function from 1. to rsc_psi*zwall/rsc_psi0
@@ -281,28 +277,21 @@
                ! building the matrix
                zcof = rfact_tke * tmask(ji,jj,jk)
-               !
-               ! lower diagonal
-               z_elem_a(ji,jj,jk) = zcof * ( avm  (ji,jj,jk  ) + avm  (ji,jj,jk-1) )   &
-                  &                      / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk  ) )
-               !
-               ! upper diagonal
-               z_elem_c(ji,jj,jk) = zcof * ( avm  (ji,jj,jk+1) + avm  (ji,jj,jk  ) )   &
-                  &                      / ( e3t_n(ji,jj,jk  ) * e3w_n(ji,jj,jk) )
-               !
-               ! diagonal
-               z_elem_b(ji,jj,jk) = 1._wp - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk)  &
-                  &                       + rdt * zdiss * tmask(ji,jj,jk) 
-               !
-               ! right hand side in en
-               en(ji,jj,jk) = en(ji,jj,jk) + rdt * zesh2 * tmask(ji,jj,jk)
+               !                                               ! lower diagonal
+               zd_lw(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk  ) + p_avm(ji,jj,jk-1) ) / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk) )
+               !                                               ! upper diagonal
+               zd_up(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk  ) ) / ( e3t_n(ji,jj,jk  ) * e3w_n(ji,jj,jk) )
+               !                                               ! diagonal
+               zdiag(ji,jj,jk) = 1._wp - zd_lw(ji,jj,jk) - zd_up(ji,jj,jk)  + rdt * zdiss * wmask(ji,jj,jk) 
+               !                                               ! right hand side in en
+               en(ji,jj,jk) = en(ji,jj,jk) + rdt * zesh2 * wmask(ji,jj,jk)
             END DO
          END DO
       END DO
       !
-      z_elem_b(:,:,jpk) = 1._wp
+      zdiag(:,:,jpk) = 1._wp
       !
       ! Set surface condition on zwall_psi (1 at the bottom)
-      zwall_psi(:,:,1) = zwall_psi(:,:,2)
-      zwall_psi(:,:,jpk) = 1.
+      zwall_psi(:,:, 1 ) = zwall_psi(:,:,2)
+      zwall_psi(:,:,jpk) = 1._wp
       !
       ! Surface boundary condition on tke
@@ -311,39 +300,36 @@
       SELECT CASE ( nn_bc_surf )
       !
-      CASE ( 0 )             ! Dirichlet case
+      CASE ( 0 )             ! Dirichlet boundary condition (set e at k=1 & 2) 
       ! First level
-      en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp)
-      en(:,:,1) = MAX(en(:,:,1), rn_emin) 
-      z_elem_a(:,:,1) = en(:,:,1)
-      z_elem_c(:,:,1) = 0._wp
-      z_elem_b(:,:,1) = 1._wp
+      en   (:,:,1) = MAX(  rn_emin , rc02r * ustar2_surf(:,:) * (1._wp + rsbc_tke1)**r2_3  )
+      zd_lw(:,:,1) = en(:,:,1)
+      zd_up(:,:,1) = 0._wp
+      zdiag(:,:,1) = 1._wp
       ! 
       ! One level below
-      en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1 * ((zhsro(:,:)+gdepw_n(:,:,2)) &
-         &               / zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp)
-      en(:,:,2) = MAX(en(:,:,2), rn_emin )
-      z_elem_a(:,:,2) = 0._wp 
-      z_elem_c(:,:,2) = 0._wp
-      z_elem_b(:,:,2) = 1._wp
-      !
-      !
-      CASE ( 1 )             ! Neumann boundary condition on d(e)/dz
+      en   (:,:,2) =  MAX(  rc02r * ustar2_surf(:,:) * (  1._wp + rsbc_tke1 * ((zhsro(:,:)+gdepw_n(:,:,2))   &
+         &                 / zhsro(:,:) )**(1.5_wp*ra_sf)  )**(2._wp/3._wp)                      , rn_emin   )
+      zd_lw(:,:,2) = 0._wp 
+      zd_up(:,:,2) = 0._wp
+      zdiag(:,:,2) = 1._wp
+      !
+      !
+      CASE ( 1 )             ! Neumann boundary condition (set d(e)/dz)
       !
       ! Dirichlet conditions at k=1
-      en(:,:,1)       = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp)
-      en(:,:,1)       = MAX(en(:,:,1), rn_emin)      
-      z_elem_a(:,:,1) = en(:,:,1)
-      z_elem_c(:,:,1) = 0._wp
-      z_elem_b(:,:,1) = 1._wp
+      en   (:,:,1) = MAX(  rc02r * ustar2_surf(:,:) * (1._wp + rsbc_tke1)**r2_3 , rn_emin  )
+      zd_lw(:,:,1) = en(:,:,1)
+      zd_up(:,:,1) = 0._wp
+      zdiag(:,:,1) = 1._wp
       !
       ! at k=2, set de/dz=Fw
       !cbr
-      z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b
-      z_elem_a(:,:,2) = 0._wp
-      zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans*gdept_n(:,:,1)/zhsro(:,:)) ))
-      zflxs(:,:)      = rsbc_tke2 * ustars2(:,:)**1.5_wp * zkar(:,:) &
-          &                       * ((zhsro(:,:)+gdept_n(:,:,1)) / zhsro(:,:) )**(1.5_wp*ra_sf)
-
-      en(:,:,2) = en(:,:,2) + zflxs(:,:)/e3w_n(:,:,2)
+      zdiag(:,:,2) = zdiag(:,:,2) +  zd_lw(:,:,2) ! Remove zd_lw from zdiag
+      zd_lw(:,:,2) = 0._wp
+      zkar (:,:)   = (rl_sf + (vkarmn-rl_sf)*(1.-EXP(-rtrans*gdept_n(:,:,1)/zhsro(:,:)) ))
+      zflxs(:,:)   = rsbc_tke2 * ustar2_surf(:,:)**1.5_wp * zkar(:,:) &
+          &                    * (  ( zhsro(:,:)+gdept_n(:,:,1) ) / zhsro(:,:)  )**(1.5_wp*ra_sf)
+!!gm why not   :                        * ( 1._wp + gdept_n(:,:,1) / zhsro(:,:) )**(1.5_wp*ra_sf)
+      en(:,:,2) = en(:,:,2) + zflxs(:,:) / e3w_n(:,:,2)
       !
       !
@@ -356,6 +342,46 @@
       !
       CASE ( 0 )             ! Dirichlet 
-         !                      ! en(ibot) = u*^2 / Co2 and mxln(ibot) = rn_lmin
+         !                      ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = rn_lmin
          !                      ! Balance between the production and the dissipation terms
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+!!gm This means that bottom and ocean w-level above have a specified "en" value.   Sure ????
+!!   With thick deep ocean level thickness, this may be quite large, no ???
+!!   in particular in ocean cavities where top stratification can be large...
+               ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point
+               ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1
+               !
+               z_en =  MAX( rc02r * ustar2_bot(ji,jj), rn_emin )
+               !
+               ! Dirichlet condition applied at: 
+               !     Bottom level (ibot)      &      Just above it (ibotm1)   
+               zd_lw(ji,jj,ibot) = 0._wp   ;   zd_lw(ji,jj,ibotm1) = 0._wp
+               zd_up(ji,jj,ibot) = 0._wp   ;   zd_up(ji,jj,ibotm1) = 0._wp
+               zdiag(ji,jj,ibot) = 1._wp   ;   zdiag(ji,jj,ibotm1) = 1._wp
+               en   (ji,jj,ibot) = z_en    ;   en   (ji,jj,ibotm1) = z_en
+            END DO
+         END DO
+         !
+         IF( ln_isfcav) THEN     ! top boundary   (ocean cavity)
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  itop   = mikt(ji,jj)       ! k   top w-point
+                  itopp1 = mikt(ji,jj) + 1   ! k+1 1st w-point below the top one
+                  !                                                ! mask at the ocean surface points
+                  z_en = MAX( rc02r * ustar2_top(ji,jj), rn_emin ) * ( 1._wp - tmask(ji,jj,1) )
+                  !
+ !!gm TO BE VERIFIED !!!
+                  ! Dirichlet condition applied at: 
+                  !     top level (itop)         &      Just below it (itopp1)   
+                  zd_lw(ji,jj,itop) = 0._wp   ;   zd_lw(ji,jj,itopp1) = 0._wp
+                  zd_up(ji,jj,itop) = 0._wp   ;   zd_up(ji,jj,itopp1) = 0._wp
+                  zdiag(ji,jj,itop) = 1._wp   ;   zdiag(ji,jj,itopp1) = 1._wp
+                  en   (ji,jj,itop) = z_en    ;   en   (ji,jj,itopp1) = z_en
+               END DO
+            END DO
+         ENDIF
+         !
+      CASE ( 1 )             ! Neumman boundary condition
+         !                      
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -363,36 +389,31 @@
                ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1
                !
+               z_en =  MAX( rc02r * ustar2_bot(ji,jj), rn_emin )
+               !
                ! Bottom level Dirichlet condition:
-               z_elem_a(ji,jj,ibot  ) = 0._wp
-               z_elem_c(ji,jj,ibot  ) = 0._wp
-               z_elem_b(ji,jj,ibot  ) = 1._wp
-               en(ji,jj,ibot  ) = MAX( rc02r * ustarb2(ji,jj), rn_emin )
-               !
-               ! Just above last level, Dirichlet condition again
-               z_elem_a(ji,jj,ibotm1) = 0._wp
-               z_elem_c(ji,jj,ibotm1) = 0._wp
-               z_elem_b(ji,jj,ibotm1) = 1._wp
-               en(ji,jj,ibotm1) = MAX( rc02r * ustarb2(ji,jj), rn_emin ) 
-            END DO
-         END DO
-         !
-      CASE ( 1 )             ! Neumman boundary condition
-         !                      
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point
-               ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1
-               !
-               ! Bottom level Dirichlet condition:
-               z_elem_a(ji,jj,ibot) = 0._wp
-               z_elem_c(ji,jj,ibot) = 0._wp
-               z_elem_b(ji,jj,ibot) = 1._wp
-               en(ji,jj,ibot) = MAX( rc02r * ustarb2(ji,jj), rn_emin )
-               !
-               ! Just above last level: Neumann condition
-               z_elem_b(ji,jj,ibotm1) = z_elem_b(ji,jj,ibotm1) + z_elem_c(ji,jj,ibotm1)   ! Remove z_elem_c from z_elem_b
-               z_elem_c(ji,jj,ibotm1) = 0._wp
-            END DO
-         END DO
+               !     Bottom level (ibot)      &      Just above it (ibotm1)   
+               !         Dirichlet            !         Neumann
+               zd_lw(ji,jj,ibot) = 0._wp   !   ! Remove zd_up from zdiag
+               zdiag(ji,jj,ibot) = 1._wp   ;   zdiag(ji,jj,ibotm1) = zdiag(ji,jj,ibotm1) + zd_up(ji,jj,ibotm1)
+               zd_up(ji,jj,ibot) = 0._wp   ;   zd_up(ji,jj,ibotm1) = 0._wp
+            END DO
+         END DO
+         IF( ln_isfcav) THEN     ! top boundary   (ocean cavity)
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  itop   = mikt(ji,jj)       ! k   top w-point
+                  itopp1 = mikt(ji,jj) + 1   ! k+1 1st w-point below the top one
+                  !                                                ! mask at the ocean surface points
+                  z_en = MAX( rc02r * ustar2_top(ji,jj), rn_emin ) * ( 1._wp - tmask(ji,jj,1) )
+                  !
+                  ! Bottom level Dirichlet condition:
+                  !     Bottom level (ibot)      &      Just above it (ibotm1)   
+                  !         Dirichlet            !         Neumann
+                  zd_lw(ji,jj,itop) = 0._wp   !   ! Remove zd_up from zdiag
+                  zdiag(ji,jj,itop) = 1._wp   ;   zdiag(ji,jj,itopp1) = zdiag(ji,jj,itopp1) + zd_up(ji,jj,itopp1)
+                  zd_up(ji,jj,itop) = 0._wp   ;   zd_up(ji,jj,itopp1) = 0._wp
+               END DO
+            END DO
+         ENDIF
          !
       END SELECT
@@ -404,5 +425,5 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1    ! vector opt.
-               z_elem_b(ji,jj,jk) = z_elem_b(ji,jj,jk) - z_elem_a(ji,jj,jk) * z_elem_c(ji,jj,jk-1) / z_elem_b(ji,jj,jk-1)
+               zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1)
             END DO
          END DO
@@ -411,5 +432,5 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1    ! vector opt.
-               z_elem_a(ji,jj,jk) = en(ji,jj,jk) - z_elem_a(ji,jj,jk) / z_elem_b(ji,jj,jk-1) * z_elem_a(ji,jj,jk-1)
+               zd_lw(ji,jj,jk) = en(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) * zd_lw(ji,jj,jk-1)
             END DO
          END DO
@@ -418,5 +439,5 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1    ! vector opt.
-               en(ji,jj,jk) = ( z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) * en(ji,jj,jk+1) ) / z_elem_b(ji,jj,jk)
+               en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk)
             END DO
          END DO
@@ -437,5 +458,5 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  psi(ji,jj,jk)  = eb(ji,jj,jk) * mxlb(ji,jj,jk)
+                  psi(ji,jj,jk)  = eb(ji,jj,jk) * hmxl_b(ji,jj,jk)
                END DO
             END DO
@@ -455,5 +476,5 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  psi(ji,jj,jk)  = SQRT( eb(ji,jj,jk) ) / ( rc0 * mxlb(ji,jj,jk) )
+                  psi(ji,jj,jk)  = SQRT( eb(ji,jj,jk) ) / ( rc0 * hmxl_b(ji,jj,jk) )
                END DO
             END DO
@@ -464,5 +485,5 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  psi(ji,jj,jk)  = rc02 * eb(ji,jj,jk) * mxlb(ji,jj,jk)**rnn 
+                  psi(ji,jj,jk)  = rc02 * eb(ji,jj,jk) * hmxl_b(ji,jj,jk)**rnn 
                END DO
             END DO
@@ -475,5 +496,5 @@
       ! 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 ).
-      ! z_elem_b : diagonal z_elem_c : upper diagonal z_elem_a : lower diagonal
+      ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal
       ! Warning : after this step, en : right hand side of the matrix
 
@@ -485,42 +506,38 @@
                zratio = psi(ji,jj,jk) / eb(ji,jj,jk) 
                !
-               ! psi3+ : stable : B=-KhN²<0 => N²>0 if rn2>0 dir = 1 (stable) otherwise dir = 0 (unstable)
-               dir = 0.5_wp + SIGN( 0.5_wp, rn2(ji,jj,jk) )
-               !
-               rpsi3 = dir * rpsi3m + ( 1._wp - dir ) * rpsi3p
+               ! psi3+ : stable : B=-KhN²<0 => N²>0 if rn2>0 zdir = 1 (stable) otherwise zdir = 0 (unstable)
+               zdir = 0.5_wp + SIGN( 0.5_wp, rn2(ji,jj,jk) )
+               !
+               rpsi3 = zdir * rpsi3m + ( 1._wp - zdir ) * rpsi3p
                !
                ! shear prod. - stratif. destruction
-               prod = rpsi1 * zratio * shear(ji,jj,jk)
+               prod = rpsi1 * zratio * p_sh2(ji,jj,jk)
                !
                ! stratif. destruction
-               buoy = rpsi3 * zratio * (- avt(ji,jj,jk) * rn2(ji,jj,jk) )
+               buoy = rpsi3 * zratio * (- p_avt(ji,jj,jk) * rn2(ji,jj,jk) )
                !
                ! shear prod. - stratif. destruction
                diss = rpsi2 * zratio * zwall(ji,jj,jk) * eps(ji,jj,jk)
                !
-               dir = 0.5_wp + SIGN( 0.5_wp, prod + buoy )   ! dir =1(=0) if shear(ji,jj,jk)+buoy >0(<0)
-               !
-               zesh2 = dir * ( prod + buoy )          + (1._wp - dir ) * prod                        ! production term
-               zdiss = dir * ( diss / psi(ji,jj,jk) ) + (1._wp - dir ) * (diss-buoy) / psi(ji,jj,jk) ! dissipation term
+               zdir = 0.5_wp + SIGN( 0.5_wp, prod + buoy )     ! zdir =1(=0) if shear(ji,jj,jk)+buoy >0(<0)
+               !
+               zesh2 = zdir * ( prod + buoy )          + (1._wp - zdir ) * prod                        ! production term
+               zdiss = zdir * ( diss / psi(ji,jj,jk) ) + (1._wp - zdir ) * (diss-buoy) / psi(ji,jj,jk) ! dissipation term
                !                                                        
                ! building the matrix
                zcof = rfact_psi * zwall_psi(ji,jj,jk) * tmask(ji,jj,jk)
-               ! lower diagonal
-               z_elem_a(ji,jj,jk) = zcof * ( avm  (ji,jj,jk  ) + avm  (ji,jj,jk-1) )   &
-                  &                      / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk  ) )
-               ! upper diagonal
-               z_elem_c(ji,jj,jk) = zcof * ( avm  (ji,jj,jk+1) + avm  (ji,jj,jk  ) )   &
-                  &                      / ( e3t_n(ji,jj,jk  ) * e3w_n(ji,jj,jk) )
-               ! diagonal
-               z_elem_b(ji,jj,jk) = 1._wp - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk)  &
-                  &                       + rdt * zdiss * tmask(ji,jj,jk)
-               !
-               ! right hand side in psi
-               psi(ji,jj,jk) = psi(ji,jj,jk) + rdt * zesh2 * tmask(ji,jj,jk)
+               !                                               ! lower diagonal
+               zd_lw(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk  ) + p_avm(ji,jj,jk-1) ) / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk) )
+               !                                               ! upper diagonal
+               zd_up(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk  ) ) / ( e3t_n(ji,jj,jk  ) * e3w_n(ji,jj,jk) )
+               !                                               ! diagonal
+               zdiag(ji,jj,jk) = 1._wp - zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) + rdt * zdiss * wmask(ji,jj,jk)
+               !                                               ! right hand side in psi
+               psi(ji,jj,jk) = psi(ji,jj,jk) + rdt * zesh2 * wmask(ji,jj,jk)
             END DO
          END DO
       END DO
       !
-      z_elem_b(:,:,jpk) = 1._wp
+      zdiag(:,:,jpk) = 1._wp
 
       ! Surface boundary condition on psi
@@ -530,45 +547,42 @@
       !
       CASE ( 0 )             ! Dirichlet boundary conditions
-      !
-      ! Surface value
-      zdep(:,:)       = zhsro(:,:) * rl_sf ! Cosmetic
-      psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
-      z_elem_a(:,:,1) = psi(:,:,1)
-      z_elem_c(:,:,1) = 0._wp
-      z_elem_b(:,:,1) = 1._wp
-      !
-      ! One level below
-      zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*gdepw_n(:,:,2)/zhsro(:,:) )))
-      zdep(:,:)       = (zhsro(:,:) + gdepw_n(:,:,2)) * zkar(:,:)
-      psi (:,:,2)     = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
-      z_elem_a(:,:,2) = 0._wp
-      z_elem_c(:,:,2) = 0._wp
-      z_elem_b(:,:,2) = 1._wp
-      ! 
-      !
+         !
+         ! Surface value
+         zdep    (:,:)   = zhsro(:,:) * rl_sf ! Cosmetic
+         psi     (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
+         zd_lw(:,:,1) = psi(:,:,1)
+         zd_up(:,:,1) = 0._wp
+         zdiag(:,:,1) = 1._wp
+         !
+         ! One level below
+         zkar    (:,:)   = (rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdepw_n(:,:,2)/zhsro(:,:) )))
+         zdep    (:,:)   = (zhsro(:,:) + gdepw_n(:,:,2)) * zkar(:,:)
+         psi     (:,:,2) = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
+         zd_lw(:,:,2) = 0._wp
+         zd_up(:,:,2) = 0._wp
+         zdiag(:,:,2) = 1._wp
+         ! 
       CASE ( 1 )             ! Neumann boundary condition on d(psi)/dz
-      !
-      ! Surface value: Dirichlet
-      zdep(:,:)       = zhsro(:,:) * rl_sf
-      psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
-      z_elem_a(:,:,1) = psi(:,:,1)
-      z_elem_c(:,:,1) = 0._wp
-      z_elem_b(:,:,1) = 1._wp
-      !
-      ! Neumann condition at k=2
-      z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b
-      z_elem_a(:,:,2) = 0._wp
-      !
-      ! Set psi vertical flux at the surface:
-      zkar(:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*gdept_n(:,:,1)/zhsro(:,:) )) ! Lengh scale slope
-      zdep(:,:) = ((zhsro(:,:) + gdept_n(:,:,1)) / zhsro(:,:))**(rmm*ra_sf)
-      zflxs(:,:) = (rnn + rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(:,:))*(1._wp + rsbc_tke1*zdep(:,:))**(2._wp*rmm/3._wp-1_wp)
-      zdep(:,:) =  rsbc_psi1 * (zwall_psi(:,:,1)*avm(:,:,1)+zwall_psi(:,:,2)*avm(:,:,2)) * &
-             & ustars2(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + gdept_n(:,:,1))**(rnn-1.)
-      zflxs(:,:) = zdep(:,:) * zflxs(:,:)
-      psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / e3w_n(:,:,2)
-
-      !   
-      !
+         !
+         ! Surface value: Dirichlet
+         zdep    (:,:)   = zhsro(:,:) * rl_sf
+         psi     (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
+         zd_lw(:,:,1) = psi(:,:,1)
+         zd_up(:,:,1) = 0._wp
+         zdiag(:,:,1) = 1._wp
+         !
+         ! Neumann condition at k=2
+         zdiag(:,:,2) = zdiag(:,:,2) +  zd_lw(:,:,2) ! Remove zd_lw from zdiag
+         zd_lw(:,:,2) = 0._wp
+         !
+         ! Set psi vertical flux at the surface:
+         zkar (:,:)   = rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdept_n(:,:,1)/zhsro(:,:) )) ! Lengh scale slope
+         zdep (:,:)   = ((zhsro(:,:) + gdept_n(:,:,1)) / zhsro(:,:))**(rmm*ra_sf)
+         zflxs(:,:)   = (rnn + rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(:,:))*(1._wp + rsbc_tke1*zdep(:,:))**(2._wp*rmm/3._wp-1_wp)
+         zdep (:,:)   = rsbc_psi1 * (zwall_psi(:,:,1)*avm(:,:,1)+zwall_psi(:,:,2)*avm(:,:,2)) * &
+            &           ustar2_surf(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + gdept_n(:,:,1))**(rnn-1.)
+         zflxs(:,:)   = zdep(:,:) * zflxs(:,:)
+         psi  (:,:,2) = psi(:,:,2) + zflxs(:,:) / e3w_n(:,:,2)
+         !
       END SELECT
 
@@ -576,9 +590,11 @@
       ! --------------------------------
       !
-      SELECT CASE ( nn_bc_bot )
-      !
+!!gm should be done for ISF (top boundary cond.)
+!!gm so, totally new staff needed      ===>>> think about that !
+!
+      SELECT CASE ( nn_bc_bot )     ! bottom boundary
       !
       CASE ( 0 )             ! Dirichlet 
-         !                      ! en(ibot) = u*^2 / Co2 and mxln(ibot) = vkarmn * rn_bfrz0
+         !                      ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = vkarmn * r_z0_bot
          !                      ! Balance between the production and the dissipation terms
          DO jj = 2, jpjm1
@@ -586,16 +602,16 @@
                ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point
                ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1
-               zdep(ji,jj) = vkarmn * rn_bfrz0
+               zdep(ji,jj) = vkarmn * r_z0_bot
                psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn
-               z_elem_a(ji,jj,ibot) = 0._wp
-               z_elem_c(ji,jj,ibot) = 0._wp
-               z_elem_b(ji,jj,ibot) = 1._wp
+               zd_lw(ji,jj,ibot) = 0._wp
+               zd_up(ji,jj,ibot) = 0._wp
+               zdiag(ji,jj,ibot) = 1._wp
                !
                ! Just above last level, Dirichlet condition again (GOTM like)
-               zdep(ji,jj) = vkarmn * ( rn_bfrz0 + e3t_n(ji,jj,ibotm1) )
+               zdep(ji,jj) = vkarmn * ( r_z0_bot + e3t_n(ji,jj,ibotm1) )
                psi (ji,jj,ibotm1) = rc0**rpp * en(ji,jj,ibot  )**rmm * zdep(ji,jj)**rnn
-               z_elem_a(ji,jj,ibotm1) = 0._wp
-               z_elem_c(ji,jj,ibotm1) = 0._wp
-               z_elem_b(ji,jj,ibotm1) = 1._wp
+               zd_lw(ji,jj,ibotm1) = 0._wp
+               zd_up(ji,jj,ibotm1) = 0._wp
+               zdiag(ji,jj,ibotm1) = 1._wp
             END DO
          END DO
@@ -609,18 +625,18 @@
                !
                ! Bottom level Dirichlet condition:
-               zdep(ji,jj) = vkarmn * rn_bfrz0
+               zdep(ji,jj) = vkarmn * r_z0_bot
                psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn
                !
-               z_elem_a(ji,jj,ibot) = 0._wp
-               z_elem_c(ji,jj,ibot) = 0._wp
-               z_elem_b(ji,jj,ibot) = 1._wp
+               zd_lw(ji,jj,ibot) = 0._wp
+               zd_up(ji,jj,ibot) = 0._wp
+               zdiag(ji,jj,ibot) = 1._wp
                !
                ! Just above last level: Neumann condition with flux injection
-               z_elem_b(ji,jj,ibotm1) = z_elem_b(ji,jj,ibotm1) + z_elem_c(ji,jj,ibotm1) ! Remove z_elem_c from z_elem_b
-               z_elem_c(ji,jj,ibotm1) = 0.
+               zdiag(ji,jj,ibotm1) = zdiag(ji,jj,ibotm1) + zd_up(ji,jj,ibotm1) ! Remove zd_up from zdiag
+               zd_up(ji,jj,ibotm1) = 0.
                !
                ! Set psi vertical flux at the bottom:
-               zdep(ji,jj) = rn_bfrz0 + 0.5_wp*e3t_n(ji,jj,ibotm1)
-               zflxb = rsbc_psi2 * ( avm(ji,jj,ibot) + avm(ji,jj,ibotm1) )   &
+               zdep(ji,jj) = r_z0_bot + 0.5_wp*e3t_n(ji,jj,ibotm1)
+               zflxb = rsbc_psi2 * ( p_avm(ji,jj,ibot) + p_avm(ji,jj,ibotm1) )   &
                   &  * (0.5_wp*(en(ji,jj,ibot)+en(ji,jj,ibotm1)))**rmm * zdep(ji,jj)**(rnn-1._wp)
                psi(ji,jj,ibotm1) = psi(ji,jj,ibotm1) + zflxb / e3w_n(ji,jj,ibotm1)
@@ -636,5 +652,5 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1    ! vector opt.
-               z_elem_b(ji,jj,jk) = z_elem_b(ji,jj,jk) - z_elem_a(ji,jj,jk) * z_elem_c(ji,jj,jk-1) / z_elem_b(ji,jj,jk-1)
+               zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1)
             END DO
          END DO
@@ -643,5 +659,5 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1    ! vector opt.
-               z_elem_a(ji,jj,jk) = psi(ji,jj,jk) - z_elem_a(ji,jj,jk) / z_elem_b(ji,jj,jk-1) * z_elem_a(ji,jj,jk-1)
+               zd_lw(ji,jj,jk) = psi(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) * zd_lw(ji,jj,jk-1)
             END DO
          END DO
@@ -650,5 +666,5 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1    ! vector opt.
-               psi(ji,jj,jk) = ( z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) * psi(ji,jj,jk+1) ) / z_elem_b(ji,jj,jk)
+               psi(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * psi(ji,jj,jk+1) ) / zdiag(ji,jj,jk)
             END DO
          END DO
@@ -703,13 +719,13 @@
       ! Limit dissipation rate under stable stratification
       ! --------------------------------------------------
-      DO jk = 1, jpkm1 ! Note that this set boundary conditions on mxln at the same time
+      DO jk = 1, jpkm1 ! Note that this set boundary conditions on hmxl_n at the same time
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1    ! vector opt.
                ! limitation
-               eps(ji,jj,jk)  = MAX( eps(ji,jj,jk), rn_epsmin )
-               mxln(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / eps(ji,jj,jk)
+               eps   (ji,jj,jk)  = MAX( eps(ji,jj,jk), rn_epsmin )
+               hmxl_n(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / eps(ji,jj,jk)
                ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated) 
                zrn2 = MAX( rn2(ji,jj,jk), rsmall )
-               IF (ln_length_lim) mxln(ji,jj,jk) = MIN(  rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), mxln(ji,jj,jk) )
+               IF( ln_length_lim )   hmxl_n(ji,jj,jk) = MIN(  rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), hmxl_n(ji,jj,jk) )
             END DO
          END DO
@@ -727,5 +743,5 @@
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   ! zcof =  l²/q²
-                  zcof = mxlb(ji,jj,jk) * mxlb(ji,jj,jk) / ( 2._wp*eb(ji,jj,jk) )
+                  zcof = hmxl_b(ji,jj,jk) * hmxl_b(ji,jj,jk) / ( 2._wp*eb(ji,jj,jk) )
                   ! Gh = -N²l²/q²
                   gh = - rn2(ji,jj,jk) * zcof
@@ -736,7 +752,7 @@
                   sm = ( rb1**(-1._wp/3._wp) + ( 18._wp*ra1*ra1 + 9._wp*ra1*ra2*(1._wp-rc2) )*sh*gh ) / (1._wp-9._wp*ra1*ra2*gh)
                   !
-                  ! Store stability function in avmu and avmv
-                  avmu(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk)
-                  avmv(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk)
+                  ! Store stability function in zstt and zstm
+                  zstt(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk)
+                  zstm(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk)
                END DO
             END DO
@@ -748,5 +764,5 @@
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   ! zcof =  l²/q²
-                  zcof = mxlb(ji,jj,jk)*mxlb(ji,jj,jk) / ( 2._wp * eb(ji,jj,jk) )
+                  zcof = hmxl_b(ji,jj,jk)*hmxl_b(ji,jj,jk) / ( 2._wp * eb(ji,jj,jk) )
                   ! Gh = -N²l²/q²
                   gh = - rn2(ji,jj,jk) * zcof
@@ -755,5 +771,5 @@
                   gh = gh * rf6
                   ! Gm =  M²l²/q² Shear number
-                  shr = shear(ji,jj,jk) / MAX( avm(ji,jj,jk), rsmall )
+                  shr = p_sh2(ji,jj,jk) / MAX( p_avm(ji,jj,jk), rsmall )
                   gm = MAX( shr * zcof , 1.e-10 )
                   gm = gm * rf6
@@ -764,7 +780,7 @@
                   sh = (rs4 - rs5*gh + rs6*gm) / rcff
                   !
-                  ! Store stability function in avmu and avmv
-                  avmu(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk)
-                  avmv(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk)
+                  ! Store stability function in zstt and zstm
+                  zstt(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk)
+                  zstm(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk)
                END DO
             END DO
@@ -776,11 +792,13 @@
       ! Lines below are useless if GOTM style Dirichlet conditions are used
 
-      avmv(:,:,1) = avmv(:,:,2)
+      zstm(:,:,1) = zstm(:,:,2)
 
       DO jj = 2, jpjm1
          DO ji = fs_2, fs_jpim1   ! vector opt.
-            avmv(ji,jj,mbkt(ji,jj)+1) = avmv(ji,jj,mbkt(ji,jj))
+            zstm(ji,jj,mbkt(ji,jj)+1) = zstm(ji,jj,mbkt(ji,jj))
          END DO
       END DO
+!!gm should be done for ISF (top boundary cond.)
+!!gm so, totally new staff needed!!gm
 
       ! Compute diffusivities/viscosities
@@ -789,44 +807,20 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zsqen         = SQRT( 2._wp * en(ji,jj,jk) ) * mxln(ji,jj,jk)
-               zav           = zsqen * avmu(ji,jj,jk)
-               avt(ji,jj,jk) = MAX( zav, avtb(jk) )*tmask(ji,jj,jk) ! apply mask for zdfmxl routine
-               zav           = zsqen * avmv(ji,jj,jk)
-               avm(ji,jj,jk) = MAX( zav, avmb(jk) ) ! Note that avm is not masked at the surface and the bottom
+               zsqen = SQRT( 2._wp * en(ji,jj,jk) ) * hmxl_n(ji,jj,jk)
+               zavt  = zsqen * zstt(ji,jj,jk)
+               zavm  = zsqen * zstm(ji,jj,jk)
+               p_avt(ji,jj,jk) = MAX( zavt, avtb(jk) ) * wmask(ji,jj,jk) ! apply mask for zdfmxl routine
+               p_avm(ji,jj,jk) = MAX( zavm, avmb(jk) )                   ! Note that avm is not masked at the surface and the bottom
             END DO
          END DO
       END DO
-      !
-      ! Lateral boundary conditions (sign unchanged)
       avt(:,:,1)  = 0._wp
-      CALL lbc_lnk( avm, 'W', 1. )   ;   CALL lbc_lnk( avt, 'W', 1. )
-
-      DO jk = 2, jpkm1            !* vertical eddy viscosity at u- and v-points
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj  ,jk) ) * umask(ji,jj,jk)
-               avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji  ,jj+1,jk) ) * vmask(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-      avmu(:,:,1) = 0._wp             ;   avmv(:,:,1) = 0._wp                 ! set surface to zero
-      CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )       ! Lateral boundary conditions
-
+      !
       IF(ln_ctl) THEN
-         CALL prt_ctl( tab3d_1=en  , clinfo1=' gls  - e: ', tab3d_2=avt, clinfo2=' t: ', ovlap=1, kdim=jpk)
-         CALL prt_ctl( tab3d_1=avmu, clinfo1=' gls  - u: ', mask1=umask,                   &
-            &          tab3d_2=avmv, clinfo2=       ' v: ', mask2=vmask, ovlap=1, kdim=jpk )
+         CALL prt_ctl( tab3d_1=en , clinfo1=' gls  - e: ', tab3d_2=avt, clinfo2=' t: ', ovlap=1, kdim=jpk)
+         CALL prt_ctl( tab3d_1=avm, clinfo1=' gls  - m: ', ovlap=1, kdim=jpk )
       ENDIF
       !
-      avt_k (:,:,:) = avt (:,:,:)
-      avm_k (:,:,:) = avm (:,:,:)
-      avmu_k(:,:,:) = avmu(:,:,:)
-      avmv_k(:,:,:) = avmv(:,:,:)
-      !
-      CALL wrk_dealloc( jpi,jpj,       zdep, zkar, zflxs, zhsro )
-      CALL wrk_dealloc( jpi,jpj,jpk,   eb, mxlb, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi )
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('zdf_gls')
-      !
+      IF( ln_timing )   CALL timing_stop('zdf_gls')
       !
    END SUBROUTINE zdf_gls
@@ -838,8 +832,7 @@
       !!                     
       !! ** Purpose :   Initialization of the vertical eddy diffivity and 
-      !!      viscosity when using a gls turbulent closure scheme
+      !!              viscosity computed using a GLS turbulent closure scheme
       !!
       !! ** Method  :   Read the namzdf_gls namelist and check the parameters
-      !!      called at the first timestep (nit000)
       !!
       !! ** input   :   Namlist namzdf_gls
@@ -848,7 +841,4 @@
       !!
       !!----------------------------------------------------------------------
-      USE dynzdf_exp
-      USE trazdf_exp
-      !
       INTEGER ::   jk    ! dummy loop indices
       INTEGER ::   ios   ! Local integer output status for namelist read
@@ -862,5 +852,5 @@
       !!----------------------------------------------------------
       !
-      IF( nn_timing == 1 )  CALL timing_start('zdf_gls_init')
+      IF( ln_timing )   CALL timing_start('zdf_gls_init')
       !
       REWIND( numnam_ref )              ! Namelist namzdf_gls in reference namelist : Vertical eddy diffivity and viscosity using gls turbulent closure scheme
@@ -875,5 +865,5 @@
       IF(lwp) THEN                     !* Control print
          WRITE(numout,*)
-         WRITE(numout,*) 'zdf_gls_init : gls turbulent closure scheme'
+         WRITE(numout,*) 'zdf_gls_init : GLS turbulent closure scheme'
          WRITE(numout,*) '~~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namzdf_gls : set gls mixing parameters'
@@ -892,17 +882,21 @@
          WRITE(numout,*) '      Type of closure                               nn_clos        = ', nn_clos
          WRITE(numout,*) '      Surface roughness (m)                         rn_hsro        = ', rn_hsro
-         WRITE(numout,*) '      Bottom roughness (m) (nambfr namelist)        rn_bfrz0       = ', rn_bfrz0
+         WRITE(numout,*)
+         WRITE(numout,*) '   Namelist namdrg_top/_bot:   used values:'
+         WRITE(numout,*) '      top    ocean cavity roughness (m)             rn_z0(_top)   = ', r_z0_top
+         WRITE(numout,*) '      Bottom seafloor     roughness (m)             rn_z0(_bot)   = ', r_z0_bot
+         WRITE(numout,*)
       ENDIF
 
-      !                                !* allocate gls arrays
+      !                                !* allocate GLS arrays
       IF( zdf_gls_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_gls_init : unable to allocate arrays' )
 
       !                                !* Check of some namelist values
-      IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_bc_surf is 0 or 1' )
-      IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_bc_surf is 0 or 1' )
-      IF( nn_z0_met < 0 .OR. nn_z0_met > 3 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_z0_met is 0, 1, 2 or 3' )
-      IF( nn_z0_met == 3 .AND. .NOT.ln_wave ) CALL ctl_stop( 'zdf_gls_init: nn_z0_met=3 requires ln_wave=T' )
-      IF( nn_stab_func  < 0 .OR. nn_stab_func  > 3 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_stab_func is 0, 1, 2 and 3' )
-      IF( nn_clos       < 0 .OR. nn_clos       > 3 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_clos is 0, 1, 2 or 3' )
+      IF( nn_bc_surf < 0   .OR. nn_bc_surf   > 1 )   CALL ctl_stop( 'zdf_gls_init: bad flag: nn_bc_surf is 0 or 1' )
+      IF( nn_bc_surf < 0   .OR. nn_bc_surf   > 1 )   CALL ctl_stop( 'zdf_gls_init: bad flag: nn_bc_surf is 0 or 1' )
+      IF( nn_z0_met  < 0   .OR. nn_z0_met    > 3 )   CALL ctl_stop( 'zdf_gls_init: bad flag: nn_z0_met is 0, 1, 2 or 3' )
+      IF( nn_z0_met == 3  .AND. .NOT.ln_wave     )   CALL ctl_stop( 'zdf_gls_init: nn_z0_met=3 requires ln_wave=T' )
+      IF( nn_stab_func < 0 .OR. nn_stab_func > 3 )   CALL ctl_stop( 'zdf_gls_init: bad flag: nn_stab_func is 0, 1, 2 and 3' )
+      IF( nn_clos      < 0 .OR. nn_clos      > 3 )   CALL ctl_stop( 'zdf_gls_init: bad flag: nn_clos is 0, 1, 2 or 3' )
 
       SELECT CASE ( nn_clos )          !* set the parameters for the chosen closure
@@ -910,5 +904,6 @@
       CASE( 0 )                              ! k-kl  (Mellor-Yamada)
          !
-         IF(lwp) WRITE(numout,*) 'The choosen closure is k-kl closed to the classical Mellor-Yamada'
+         IF(lwp) WRITE(numout,*) '   ==>>   k-kl closure chosen (i.e. closed to the classical Mellor-Yamada)'
+         IF(lwp) WRITE(numout,*)
          rpp     = 0._wp
          rmm     = 1._wp
@@ -928,5 +923,6 @@
       CASE( 1 )                              ! k-eps
          !
-         IF(lwp) WRITE(numout,*) 'The choosen closure is k-eps'
+         IF(lwp) WRITE(numout,*) '   ==>>   k-eps closure chosen'
+         IF(lwp) WRITE(numout,*)
          rpp     =  3._wp
          rmm     =  1.5_wp
@@ -946,5 +942,6 @@
       CASE( 2 )                              ! k-omega
          !
-         IF(lwp) WRITE(numout,*) 'The choosen closure is k-omega'
+         IF(lwp) WRITE(numout,*) '   ==>>   k-omega closure chosen'
+         IF(lwp) WRITE(numout,*)
          rpp     = -1._wp
          rmm     =  0.5_wp
@@ -964,5 +961,6 @@
       CASE( 3 )                              ! generic
          !
-         IF(lwp) WRITE(numout,*) 'The choosen closure is generic'
+         IF(lwp) WRITE(numout,*) '   ==>>   generic closure chosen'
+         IF(lwp) WRITE(numout,*)
          rpp     = 2._wp
          rmm     = 1._wp
@@ -987,5 +985,5 @@
       CASE ( 0 )                             ! Galperin stability functions
          !
-         IF(lwp) WRITE(numout,*) 'Stability functions from Galperin'
+         IF(lwp) WRITE(numout,*) '   ==>>   Stability functions from Galperin'
          rc2     =  0._wp
          rc3     =  0._wp
@@ -999,5 +997,5 @@
       CASE ( 1 )                             ! Kantha-Clayson stability functions
          !
-         IF(lwp) WRITE(numout,*) 'Stability functions from Kantha-Clayson'
+         IF(lwp) WRITE(numout,*) '   ==>>   Stability functions from Kantha-Clayson'
          rc2     =  0.7_wp
          rc3     =  0.2_wp
@@ -1011,5 +1009,5 @@
       CASE ( 2 )                             ! Canuto A stability functions
          !
-         IF(lwp) WRITE(numout,*) 'Stability functions from Canuto A'
+         IF(lwp) WRITE(numout,*) '   ==>>   Stability functions from Canuto A'
          rs0 = 1.5_wp * rl1 * rl5*rl5
          rs1 = -rl4*(rl6+rl7) + 2._wp*rl4*rl5*(rl1-(1._wp/3._wp)*rl2-rl3) + 1.5_wp*rl1*rl5*rl8
@@ -1035,5 +1033,5 @@
       CASE ( 3 )                             ! Canuto B stability functions
          !
-         IF(lwp) WRITE(numout,*) 'Stability functions from Canuto B'
+         IF(lwp) WRITE(numout,*) '   ==>>   Stability functions from Canuto B'
          rs0 = 1.5_wp * rm1 * rm5*rm5
          rs1 = -rm4 * (rm6+rm7) + 2._wp * rm4*rm5*(rm1-(1._wp/3._wp)*rm2-rm3) + 1.5_wp * rm1*rm5*rm8
@@ -1079,10 +1077,9 @@
          rl_sf = vkarmn
       ELSE
-         rl_sf = rc0 * SQRT(rc0/rcm_sf) * SQRT( ( (1._wp + 4._wp*rmm + 8._wp*rmm**2_wp)*rsc_tke          &
-                 &                                       + 12._wp * rsc_psi0*rpsi2 - (1._wp + 4._wp*rmm) &
-                 &                                                *SQRT(rsc_tke*(rsc_tke                 &
-                 &                                                   + 24._wp*rsc_psi0*rpsi2)) )         &
-                 &                                         /(12._wp*rnn**2.)                             &
-                 &                                       )
+         rl_sf = rc0 * SQRT(rc0/rcm_sf) * SQRT( ( (1._wp + 4._wp*rmm + 8._wp*rmm**2_wp) * rsc_tke        &
+            &                                            + 12._wp*rsc_psi0*rpsi2 - (1._wp + 4._wp*rmm)   &
+            &                                                     *SQRT(rsc_tke*(rsc_tke                 &
+            &                                                        + 24._wp*rsc_psi0*rpsi2)) )         &
+            &                                              /(12._wp*rnn**2.)                             )
       ENDIF
 
@@ -1090,23 +1087,21 @@
       IF(lwp) THEN                     !* Control print
          WRITE(numout,*)
-         WRITE(numout,*) 'Limit values'
-         WRITE(numout,*) '~~~~~~~~~~~~'
-         WRITE(numout,*) 'Parameter  m = ',rmm
-         WRITE(numout,*) 'Parameter  n = ',rnn
-         WRITE(numout,*) 'Parameter  p = ',rpp
-         WRITE(numout,*) 'rpsi1   = ',rpsi1
-         WRITE(numout,*) 'rpsi2   = ',rpsi2
-         WRITE(numout,*) 'rpsi3m  = ',rpsi3m
-         WRITE(numout,*) 'rpsi3p  = ',rpsi3p
-         WRITE(numout,*) 'rsc_tke = ',rsc_tke
-         WRITE(numout,*) 'rsc_psi = ',rsc_psi
-         WRITE(numout,*) 'rsc_psi0 = ',rsc_psi0
-         WRITE(numout,*) 'rc0     = ',rc0
+         WRITE(numout,*) '   Limit values :'
+         WRITE(numout,*) '      Parameter  m = ', rmm
+         WRITE(numout,*) '      Parameter  n = ', rnn
+         WRITE(numout,*) '      Parameter  p = ', rpp
+         WRITE(numout,*) '      rpsi1    = ', rpsi1
+         WRITE(numout,*) '      rpsi2    = ', rpsi2
+         WRITE(numout,*) '      rpsi3m   = ', rpsi3m
+         WRITE(numout,*) '      rpsi3p   = ', rpsi3p
+         WRITE(numout,*) '      rsc_tke  = ', rsc_tke
+         WRITE(numout,*) '      rsc_psi  = ', rsc_psi
+         WRITE(numout,*) '      rsc_psi0 = ', rsc_psi0
+         WRITE(numout,*) '      rc0      = ', rc0
          WRITE(numout,*)
-         WRITE(numout,*) 'Shear free turbulence parameters:'
-         WRITE(numout,*) 'rcm_sf  = ',rcm_sf
-         WRITE(numout,*) 'ra_sf   = ',ra_sf
-         WRITE(numout,*) 'rl_sf   = ',rl_sf
-         WRITE(numout,*)
+         WRITE(numout,*) '   Shear free turbulence parameters:'
+         WRITE(numout,*) '      rcm_sf   = ', rcm_sf
+         WRITE(numout,*) '      ra_sf    = ', ra_sf
+         WRITE(numout,*) '      rl_sf    = ', rl_sf
       ENDIF
 
@@ -1123,22 +1118,16 @@
       rsbc_psi1 = -0.5_wp * rdt * rc0**(rpp-2._wp*rmm) / rsc_psi
       rsbc_psi2 = -0.5_wp * rdt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking 
-
+      !
       rfact_tke = -0.5_wp / rsc_tke * rdt                                ! Cst used for the Diffusion term of tke
       rfact_psi = -0.5_wp / rsc_psi * rdt                                ! Cst used for the Diffusion term of tke
-
+      !
       !                                !* Wall proximity function
-      zwall (:,:,:) = 1._wp * tmask(:,:,:)
-
-      !                                !* set vertical eddy coef. to the background value
-      DO jk = 1, jpk
-         avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
-         avm (:,:,jk) = avmb(jk) * tmask(:,:,jk)
-         avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)
-         avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)
-      END DO
-      !                              
-      CALL gls_rst( nit000, 'READ' )   !* read or initialize all required files
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('zdf_gls_init')
+!!gm tmask or wmask ????
+      zwall(:,:,:) = 1._wp * tmask(:,:,:)
+
+      !                                !* read or initialize all required files  
+      CALL gls_rst( nit000, 'READ' )      ! (en, avt_k, avm_k, hmxl_n)
+      !
+      IF( ln_timing )   CALL timing_stop('zdf_gls_init')
       !
    END SUBROUTINE zdf_gls_init
@@ -1155,9 +1144,9 @@
       !!                set to rn_emin or recomputed (nn_igls/=0)
       !!----------------------------------------------------------------------
-      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
-      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
+      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
+      INTEGER ::   id1, id2, id3, id4
       INTEGER ::   ji, jj, ikbu, ikbv
       REAL(wp)::   cbx, cby
@@ -1167,32 +1156,27 @@
          !                                   ! ---------------
          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, 'mxln' , ldstop = .FALSE. )
+            id1 = iom_varid( numror, 'en'    , ldstop = .FALSE. )
+            id2 = iom_varid( numror, 'avt_k' , ldstop = .FALSE. )
+            id3 = iom_varid( numror, 'avm_k' , ldstop = .FALSE. )
+            id4 = iom_varid( numror, 'hmxl_n', ldstop = .FALSE. )
             !
-            IF( MIN( id1, id2, id3, id4, id5, id6 ) > 0 ) THEN        ! all required arrays exist
+            IF( MIN( id1, id2, id3, id4 ) > 0 ) THEN        ! all required arrays exist
                CALL iom_get( numror, jpdom_autoglo, 'en'    , en     )
-               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, 'mxln'  , mxln   )
+               CALL iom_get( numror, jpdom_autoglo, 'avt_k' , avt_k  )
+               CALL iom_get( numror, jpdom_autoglo, 'avm_k' , avm_k  )
+               CALL iom_get( numror, jpdom_autoglo, 'hmxl_n', hmxl_n )
             ELSE                        
-               IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without gls scheme, en and mxln computed by iterative loop'
-               en  (:,:,:) = rn_emin
-               mxln(:,:,:) = 0.05        
-               avt_k (:,:,:) = avt (:,:,:)
-               avm_k (:,:,:) = avm (:,:,:)
-               avmu_k(:,:,:) = avmu(:,:,:)
-               avmv_k(:,:,:) = avmv(:,:,:)
-               DO jit = nit000 + 1, nit000 + 10   ;   CALL zdf_gls( jit )   ;   END DO
+               IF(lwp) WRITE(numout,*)
+               IF(lwp) WRITE(numout,*) '   ==>>   previous run without GLS scheme, set en and hmxl_n to background values'
+               en    (:,:,:) = rn_emin
+               hmxl_n(:,:,:) = 0.05_wp
+               ! avt_k, avm_k already set to the background value in zdf_phy_init
             ENDIF
          ELSE                                   !* Start from rest
-            IF(lwp) WRITE(numout,*) ' ===>>>> : Initialisation of en and mxln by background values'
-            en  (:,:,:) = rn_emin
-            mxln(:,:,:) = 0.05       
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) '   ==>>   start from rest, set en and hmxl_n by background values'
+            en    (:,:,:) = rn_emin
+            hmxl_n(:,:,:) = 0.05_wp
+            ! avt_k, avm_k already set to the background value in zdf_phy_init
          ENDIF
          !
@@ -1200,33 +1184,12 @@
          !                                   ! -------------------
          IF(lwp) WRITE(numout,*) '---- gls-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, 'mxln' , mxln   )
+         CALL iom_rstput( kt, nitrst, numrow, 'en'    , en     ) 
+         CALL iom_rstput( kt, nitrst, numrow, 'avt_k' , avt_k  )
+         CALL iom_rstput( kt, nitrst, numrow, 'avm_k' , avm_k  )
+         CALL iom_rstput( kt, nitrst, numrow, 'hmxl_n', hmxl_n )
          !
       ENDIF
       !
    END SUBROUTINE gls_rst
-
-#else
-   !!----------------------------------------------------------------------
-   !!   Dummy module :                                        NO TKE scheme
-   !!----------------------------------------------------------------------
-   LOGICAL, PUBLIC, PARAMETER ::   lk_zdfgls = .FALSE.   !: TKE flag
-CONTAINS
-   SUBROUTINE zdf_gls_init           ! Empty routine
-      WRITE(*,*) 'zdf_gls_init: You should not have seen this print! error?'
-   END SUBROUTINE zdf_gls_init
-   SUBROUTINE zdf_gls( kt )          ! Empty routine
-      WRITE(*,*) 'zdf_gls: You should not have seen this print! error?', kt
-   END SUBROUTINE zdf_gls
-   SUBROUTINE gls_rst( kt, cdrw )          ! Empty routine
-      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
-      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
-      WRITE(*,*) 'gls_rst: You should not have seen this print! error?', kt, cdrw
-   END SUBROUTINE gls_rst
-#endif
 
    !!======================================================================
