Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90	(revision 7953)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90	(revision 7990)
@@ -5,6 +5,7 @@
    !!                 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 
    !!----------------------------------------------------------------------
 
@@ -36,10 +37,10 @@
    PRIVATE
 
-   PUBLIC   zdf_gls        ! routine called in step module
-   PUBLIC   zdf_gls_init   ! routine called in zdfphy module
-   PUBLIC   gls_rst        ! routine called in zdfphy module
+   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(:,:,:) ::   hmxn    !: 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
@@ -113,5 +114,5 @@
       !!                ***  FUNCTION zdf_gls_alloc  ***
       !!----------------------------------------------------------------------
-      ALLOCATE( mxln(jpi,jpj,jpk), zwall(jpi,jpj,jpk) ,     &
+      ALLOCATE( hmxn(jpi,jpj,jpk), zwall(jpi,jpj,jpk) ,     &
          &      ustars2(jpi,jpj) , ustarb2(jpi,jpj)   , STAT= zdf_gls_alloc )
          !
@@ -148,21 +149,20 @@
       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), POINTER, DIMENSION(:,:,:) ::   z3du, z3dv  ! u- and v-shears
       !!--------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )  CALL timing_start('zdf_gls')
+      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  )
-      
+      CALL wrk_alloc( jpi,jpj,jpk,   eb, mxlb, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi )
+      CALL wrk_alloc( jpi,jpj,jpk,   z3du, z3dv )
+
       ! 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
+      ! restore before values computed GLS alone
+      avt(:,:,:) = avt_k(:,:,:)
+      avm(:,:,:) = avm_k(:,:,:)
 
       ! Compute surface and bottom friction at T-points
@@ -183,39 +183,37 @@
       END DO    
 
-      ! Set surface roughness length
-      SELECT CASE ( nn_z0_met )
-      !
-      CASE ( 0 )             ! Constant roughness          
+      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)
       CASE ( 2 )             ! Roughness formulae according to Rascle et al., Ocean Modelling (2008)
+!!gm         zcof = 2._wp * 0.6_wp / 28._wp
+!!gm         zdep(:,:)  = 30._wp * TANH(  zcof/ SQRT( MAX(ustars2(:,:),rsmall) )  )      ! Wave age (eq. 10)
          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)
       CASE ( 3 )             ! Roughness given by the wave model (coupled or read in file)
+!!gm  BUG missing a multiplicative coefficient....
          zhsro(:,:) = hsw(:,:)
       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 shear and dissipation rate  ==!
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1
+               z3du(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk  ) + avm(ji+1,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) )
+               z3dv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk  ) + avm(ji,jj+1,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)) / hmxn(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(:,:,:)
+      mxlb(:,:,:) = hmxn(:,:,:)
 
       IF( nn_clos == 0 ) THEN    ! Mellor-Yamada
@@ -223,5 +221,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   = hmxn(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 ) )
@@ -247,9 +245,9 @@
       DO jk = 2, jpkm1
          DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
+            DO ji = 2, jpim1
                !
                ! 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) )
+               shear(ji,jj,jk) =  ( z3du(ji-1,jj,jk) + z3du(ji,jj,jk) ) / MAX( 1.e0 , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   &
+                  &             + ( z3dv(ji,jj-1,jk) + z3dv(ji,jj,jk) ) / MAX( 1.e0 , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )
                !
                ! stratif. destruction
@@ -278,19 +276,12 @@
                ! 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
+               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 * 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
@@ -300,6 +291,6 @@
       !
       ! 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
@@ -353,5 +344,5 @@
       !
       CASE ( 0 )             ! Dirichlet 
-         !                      ! en(ibot) = u*^2 / Co2 and mxln(ibot) = rn_lmin
+         !                      ! en(ibot) = u*^2 / Co2 and hmxn(ibot) = rn_lmin
          !                      ! Balance between the production and the dissipation terms
          DO jj = 2, jpjm1
@@ -503,16 +494,12 @@
                ! 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
+               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 * 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
@@ -565,5 +552,4 @@
       zflxs(:,:) = zdep(:,:) * zflxs(:,:)
       psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / e3w_n(:,:,2)
-
       !   
       !
@@ -577,5 +563,5 @@
       !
       CASE ( 0 )             ! Dirichlet 
-         !                      ! en(ibot) = u*^2 / Co2 and mxln(ibot) = vkarmn * rn_bfrz0
+         !                      ! en(ibot) = u*^2 / Co2 and hmxn(ibot) = vkarmn * rn_bfrz0
          !                      ! Balance between the production and the dissipation terms
          DO jj = 2, jpjm1
@@ -700,13 +686,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 hmxn 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)
+               hmxn(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) hmxn(ji,jj,jk) = MIN(  rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), hmxn(ji,jj,jk) )
             END DO
          END DO
@@ -733,7 +719,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 z3du and z3dv
+                  z3du(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk)
+                  z3dv(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk)
                END DO
             END DO
@@ -761,7 +747,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 z3du and z3dv
+                  z3du(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk)
+                  z3dv(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk)
                END DO
             END DO
@@ -773,9 +759,9 @@
       ! Lines below are useless if GOTM style Dirichlet conditions are used
 
-      avmv(:,:,1) = avmv(:,:,2)
+      z3dv(:,:,1) = z3dv(:,:,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))
+            z3dv(ji,jj,mbkt(ji,jj)+1) = z3dv(ji,jj,mbkt(ji,jj))
          END DO
       END DO
@@ -786,44 +772,31 @@
          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) ) * hmxn(ji,jj,jk)
+               zav           = zsqen * z3du(ji,jj,jk)
+               avt(ji,jj,jk) = MAX( zav, avtb(jk) ) * wmask(ji,jj,jk) ! apply mask for zdfmxl routine
+               zav           = zsqen * z3dv(ji,jj,jk)
+               avm(ji,jj,jk) = MAX( zav, 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
+!!gm I'm sure this is needed to compute the shear term at next time-step
       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(:,:,:)
+      avt_k(:,:,:) = avt(:,:,:)      !==  store avt, avm values computed by GLS only  ==!
+      avm_k(:,:,:) = avm(:,:,:)
+      !
+      IF( lrst_oce )   CALL gls_rst( kt, 'WRITE' )     ! write the GLS info in the restart file
       !
       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')
-      !
+      CALL wrk_dealloc( jpi,jpj,jpk,   z3du, z3dv )
+      !
+      IF( nn_timing == 1 )   CALL timing_stop('zdf_gls')
       !
    END SUBROUTINE zdf_gls
@@ -835,8 +808,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
@@ -892,5 +864,5 @@
       ENDIF
 
-      !                                !* allocate gls arrays
+      !                                !* allocate GLS arrays
       IF( zdf_gls_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_gls_init : unable to allocate arrays' )
 
@@ -1120,20 +1092,13 @@
       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
+      !                                !* read or initialize all required files  
+      CALL gls_rst( nit000, 'READ' )      ! (en, avt_k, avm_k, hmxn)
       !
       IF( nn_timing == 1 )  CALL timing_stop('zdf_gls_init')
@@ -1152,9 +1117,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
@@ -1165,31 +1130,29 @@
          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. )
+            id2 = iom_varid( numror, 'avt_k', ldstop = .FALSE. )
+            id3 = iom_varid( numror, 'avm_k', ldstop = .FALSE. )
+            id4 = iom_varid( numror, 'hmxn' , ldstop = .FALSE. )
             !
-            IF( MIN( id1, id2, id3, id4, id5, id6 ) > 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   )
+            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_k', avt_k )
+               CALL iom_get( numror, jpdom_autoglo, 'avm_k', avm_k )
+               CALL iom_get( numror, jpdom_autoglo, 'hmxn' , hmxn  )
             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(:,:,:)
+               IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without gls scheme, en and hmxn computed by iterative loop'
+               en   (:,:,:) = rn_emin
+               hmxn (:,:,:) = 0.05_wp
+               avt_k(:,:,:) = avt(:,:,:)
+               avm_k(:,:,:) = avm(:,:,:)
                DO jit = nit000 + 1, nit000 + 10   ;   CALL zdf_gls( jit )   ;   END DO
             ENDIF
          ELSE                                   !* Start from rest
-            IF(lwp) WRITE(numout,*) ' ===>>>> : Initialisation of en and mxln by background values'
+            IF(lwp) WRITE(numout,*) ' ===>>>> : Initialisation of en and hmxn by background values'
             en  (:,:,:) = rn_emin
-            mxln(:,:,:) = 0.05       
+            hmxn(:,:,:) = 0.05_wp
+            DO jk = 1, jpk
+               avt_k(:,:,jk) = avtb(jk) * wmask(:,:,jk)
+               avm_k(:,:,jk) = avmb(jk) * wmask(:,:,jk)
+            END DO
          ENDIF
          !
@@ -1197,10 +1160,8 @@
          !                                   ! -------------------
          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, 'hmxn' , hmxn  )
          !
       ENDIF
