Index: branches/devmercator2010_1/NEMO/OPA_SRC/DYN/dynzdf.F90
===================================================================
--- branches/devmercator2010_1/NEMO/OPA_SRC/DYN/dynzdf.F90	(revision 2131)
+++ branches/devmercator2010_1/NEMO/OPA_SRC/DYN/dynzdf.F90	(revision 2132)
@@ -107,4 +107,5 @@
       USE zdftke_old
       USE zdftke
+      USE zdfgls
       USE zdfkpp
       !!----------------------------------------------------------------------
@@ -116,7 +117,7 @@
 
       ! Force implicit schemes
-      IF( lk_zdftke_old .OR. lk_zdftke .OR. lk_zdfkpp )   nzdf = 1   ! TKE or KPP physics
-      IF( ln_dynldf_iso                               )   nzdf = 1   ! iso-neutral lateral physics
-      IF( ln_dynldf_hor .AND. ln_sco                  )   nzdf = 1   ! horizontal lateral physics in s-coordinate
+      IF( lk_zdftke_old .OR. lk_zdftke .OR. lk_zdfgls .OR. lk_zdfkpp )   nzdf = 1   ! TKE, GLS or KPP physics
+      IF( ln_dynldf_iso                                              )   nzdf = 1   ! iso-neutral lateral physics
+      IF( ln_dynldf_hor .AND. ln_sco                                 )   nzdf = 1   ! horizontal lateral physics in s-coordinate
 
       IF( lk_esopa )    nzdf = -1                   ! Esopa key: All schemes used
Index: branches/devmercator2010_1/NEMO/OPA_SRC/DYN/dynzdf_imp.F90
===================================================================
--- branches/devmercator2010_1/NEMO/OPA_SRC/DYN/dynzdf_imp.F90	(revision 2131)
+++ branches/devmercator2010_1/NEMO/OPA_SRC/DYN/dynzdf_imp.F90	(revision 2132)
@@ -24,4 +24,8 @@
    USE phycst          ! physical constants
    USE in_out_manager  ! I/O manager
+#if defined key_zdfgls
+   USE zdfbfr, ONLY : bfrua, bfrva, wbotu, wbotv ! bottom stresses
+   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
+#endif
 
    IMPLICIT NONE
@@ -75,4 +79,8 @@
       REAL(wp) ::   zzwi                               ! temporary scalars
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwi        ! temporary workspace arrays
+#if defined key_zdfgls
+      INTEGER :: ikbu, ikbv, ikbum1, ikbvm1
+      REAL(wp) :: zcbcu, zcbcv
+#endif
       !!----------------------------------------------------------------------
 
@@ -173,4 +181,16 @@
          END DO
       END DO
+
+#if defined key_zdfgls
+      ! Save bottom stress for next time step
+      DO jj = 2, jpjm1
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            ikbu   = MIN( mbathy(ji+1,jj  ), mbathy(ji,jj) )
+            ikbum1 = MAX( ikbu-1, 1 )
+            wbotu(ji,jj) = bfrua(ji,jj) * ua(ji,jj,ikbum1) * umask(ji,jj,ikbum1)
+         END DO
+      END DO
+      CALL lbc_lnk( wbotu(:,:), 'U', -1. )
+#endif
 
       ! Normalization to obtain the general momentum trend ua
@@ -272,4 +292,16 @@
       END DO
 
+#if defined key_zdfgls
+      ! Save bottom stress for next time step
+      DO jj = 2, jpjm1
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            ikbv   = MIN( mbathy(ji,jj+1), mbathy(ji,jj) )
+            ikbvm1 = MAX( ikbv-1, 1 )
+            wbotv(ji,jj) = bfrva(ji,jj) * va(ji,jj,ikbvm1) * vmask(ji,jj,ikbvm1)
+         END DO
+      END DO
+      CALL lbc_lnk( wbotv(:,:), 'V', -1. )
+#endif
+
       ! Normalization to obtain the general momentum trend va
       DO jk = 1, jpkm1
Index: branches/devmercator2010_1/NEMO/OPA_SRC/IOM/restart.F90
===================================================================
--- branches/devmercator2010_1/NEMO/OPA_SRC/IOM/restart.F90	(revision 2131)
+++ branches/devmercator2010_1/NEMO/OPA_SRC/IOM/restart.F90	(revision 2132)
@@ -26,4 +26,8 @@
    USE zdfmxl          ! mixed layer depth
    USE trdmld_oce      ! ocean active mixed layer tracers trends variables
+#if defined key_zdfgls
+   USE zdfbfr, ONLY : wbotu, wbotv ! bottom stresses
+   USE zdf_oce
+#endif
 
    IMPLICIT NONE
@@ -136,4 +140,10 @@
 #endif
 
+#if defined key_zdfgls
+      ! Save bottom stresses
+      CALL iom_rstput( kt, nitrst, numrow, 'wbotu' , wbotu )
+      CALL iom_rstput( kt, nitrst, numrow, 'wbotv' , wbotv )
+#endif
+
       IF( kt == nitrst ) THEN
          CALL iom_close( numrow )     ! close the restart file (only at last time step)
Index: branches/devmercator2010_1/NEMO/OPA_SRC/TRA/trazdf.F90
===================================================================
--- branches/devmercator2010_1/NEMO/OPA_SRC/TRA/trazdf.F90	(revision 2131)
+++ branches/devmercator2010_1/NEMO/OPA_SRC/TRA/trazdf.F90	(revision 2132)
@@ -125,4 +125,5 @@
       USE zdftke_old
       USE zdftke
+      USE zdfgls
       USE zdfkpp
       !!----------------------------------------------------------------------
@@ -139,7 +140,7 @@
 
       ! Force implicit schemes
-      IF( lk_zdftke_old .OR. lk_zdftke .OR. lk_zdfkpp )   nzdf = 1      ! TKE or KPP physics
-      IF( ln_traldf_iso                               )   nzdf = 1      ! iso-neutral lateral physics
-      IF( ln_traldf_hor .AND. ln_sco                  )   nzdf = 1      ! horizontal lateral physics in s-coordinate
+      IF( lk_zdftke_old .OR. lk_zdftke .OR. lk_zdfgls .OR. lk_zdfkpp )   nzdf = 1      ! TKE, GLS or KPP physics
+      IF( ln_traldf_iso                                              )   nzdf = 1      ! iso-neutral lateral physics
+      IF( ln_traldf_hor .AND. ln_sco                                 )   nzdf = 1      ! horizontal lateral physics in s-coordinate
 
       IF( ln_zdfexp .AND. nzdf == 1 )   THEN
Index: branches/devmercator2010_1/NEMO/OPA_SRC/ZDF/zdfbfr.F90
===================================================================
--- branches/devmercator2010_1/NEMO/OPA_SRC/ZDF/zdfbfr.F90	(revision 2131)
+++ branches/devmercator2010_1/NEMO/OPA_SRC/ZDF/zdfbfr.F90	(revision 2132)
@@ -29,4 +29,8 @@
    
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   bfrua , bfrva   !: Bottom friction coefficients set in zdfbfr
+#if defined key_zdfgls
+   REAL(wp), PUBLIC                     ::   rn_hbro =  0.003_wp  ! Bottom roughness (m)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   wbotu, wbotv    !  Bottom stresses
+#endif
 
    !                                    !!* Namelist nambfr: bottom friction namelist *
Index: branches/devmercator2010_1/NEMO/OPA_SRC/ZDF/zdfgls.F90
===================================================================
--- branches/devmercator2010_1/NEMO/OPA_SRC/ZDF/zdfgls.F90	(revision 2132)
+++ branches/devmercator2010_1/NEMO/OPA_SRC/ZDF/zdfgls.F90	(revision 2132)
@@ -0,0 +1,1342 @@
+MODULE zdfgls
+   !!======================================================================
+   !!                       ***  MODULE  zdfgls  ***
+   !! Ocean physics:  vertical mixing coefficient computed from the gls 
+   !!                 turbulent closure parameterization
+   !!======================================================================
+   !! History :   3.0  !  2009-09  (G. Reffray) : Original code
+   !!----------------------------------------------------------------------
+#if defined key_zdfgls   ||   defined key_esopa
+   !!----------------------------------------------------------------------
+   !!   'key_zdfgls'                 Generic Length Scale vertical physics
+   !!----------------------------------------------------------------------
+   !!----------------------------------------------------------------------
+   !!   zdf_gls      : update momentum and tracer Kz from a gls scheme
+   !!   zdf_gls_init : initialization, namelist read, and parameters control
+   !!   gls_rst      : read/write gls restart in ocean restart file
+   !!----------------------------------------------------------------------
+   USE oce            ! ocean dynamics and active tracers 
+   USE dom_oce        ! ocean space and time domain
+   USE domvvl         ! ocean space and time domain : variable volume layer
+   USE zdf_oce        ! ocean vertical physics
+   USE sbc_oce        ! surface boundary condition: ocean
+   USE phycst         ! physical constants
+   USE zdfmxl         ! mixed layer
+   USE restart        ! only for lrst_oce
+   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
+   USE prtctl         ! Print control
+   USE in_out_manager ! I/O manager
+   USE iom            ! I/O manager library
+   USE zdfbfr, ONLY : rn_hbro, wbotu, wbotv ! bottom roughness and bottom stresses
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   zdf_gls    ! routine called in step module
+   PUBLIC   gls_rst    ! routine called in step module
+
+   LOGICAL , PUBLIC, PARAMETER              ::   lk_zdfgls = .TRUE.  !: TKE vertical mixing flag
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   en                  !: now turbulent kinetic energy
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   mxln                !: now mixing length
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   zwall              !: wall function
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)     ::   ustars2            !: Squared surface velocity scale at T-points
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)     ::   ustarb2            !: Squared bottom  velocity scale at T-points
+
+   !                                         !!! ** Namelist  namzdf_gls  **
+   LOGICAL  ::   ln_crban   = .FALSE.         ! =T use Craig and Banner scheme
+   LOGICAL  ::   ln_length_lim = .FALSE.      ! use limit on the dissipation rate understable stratification (Galperin et al., 1988)
+   LOGICAL  ::   ln_sigpsi = .FALSE.          ! Activate Burchard (2003) modification for k-eps closure AND wave breaking mixing
+   REAL(wp) ::   rn_epsmin  = 1.e-12_wp       ! minimum value of dissipation (m2/s3)
+   REAL(wp) ::   rn_emin    = 1.e-6_wp        ! minimum value of TKE (m2/s2)
+   INTEGER  ::   nn_tkebc_surf = 0            ! TKE surface boundary condition (=0/1)
+   INTEGER  ::   nn_tkebc_bot  = 0            ! TKE bottom boundary condition (=0/1)
+   INTEGER  ::   nn_psibc_surf = 0            ! PSI surface boundary condition (=0/1)
+   INTEGER  ::   nn_psibc_bot  = 0            ! PSI bottom boundary condition (=0/1)
+   INTEGER  ::   nn_stab_func  = 0            ! stability functions G88, KC or Canuto (=0/1/2)
+   INTEGER  ::   nn_clos    = 0               ! closure 0/1/2/3 MY82/k-eps/k-w/gen
+   REAL(wp) ::   clim_galp  = 0.53_wp         ! Holt 2008 value for k-eps: 0.267
+   REAL(wp) ::   hsro       = 0.003_wp        ! Minimum surface roughness
+   REAL(wp) ::   rn_charn   = 2.e+5_wp        ! Charnock constant for surface breaking waves mixing : 1400. (standard) or 2.e5 (Stacey value)
+   REAL(wp) ::   rn_crban   = 100._wp         ! Craig and Banner constant for surface breaking waves mixing
+   REAL(wp) ::   rn_cm_sf   = 0.73_wp         ! Shear free turbulence parameters
+   REAL(wp) ::   rn_a_sf    = -2.0_wp         ! Must be negative -2 < rn_a_sf < -1 
+   REAL(wp) ::   rn_l_sf    =  0.2_wp         ! 0 <rn_l_sf<vkarmn    
+   REAL(wp) ::   rn_ghmin   = -0.28
+   REAL(wp) ::   rn_gh0     = 0.0329
+   REAL(wp) ::   rn_ghcri   = 0.03
+   REAL(wp) ::   rn_a1      =  0.92_wp
+   REAL(wp) ::   rn_a2      =  0.74_wp
+   REAL(wp) ::   rn_b1      = 16.60_wp
+   REAL(wp) ::   rn_b2      = 10.10_wp         
+   REAL(wp) ::   rn_e1      =  1.80_wp         
+   REAL(wp) ::   rn_e2      =  1.33_wp         
+   REAL(wp) ::   rn_l1      =  0.107_wp
+   REAL(wp) ::   rn_l2      =  0.0032_wp
+   REAL(wp) ::   rn_l3      =  0.0864_wp
+   REAL(wp) ::   rn_l4      =  0.12_wp
+   REAL(wp) ::   rn_l5      = 11.9_wp
+   REAL(wp) ::   rn_l6      =  0.4_wp
+   REAL(wp) ::   rn_l7      =  0.0_wp
+   REAL(wp) ::   rn_l8      =  0.48_wp
+   REAL(wp) ::   rn_m1      =  0.127_wp
+   REAL(wp) ::   rn_m2      =  0.00336_wp
+   REAL(wp) ::   rn_m3      =  0.0906_wp
+   REAL(wp) ::   rn_m4      =  0.101_wp
+   REAL(wp) ::   rn_m5      = 11.2_wp
+   REAL(wp) ::   rn_m6      =  0.4_wp
+   REAL(wp) ::   rn_m7      =  0.0_wp
+   REAL(wp) ::   rn_m8      =  0.318_wp
+   REAL(wp) ::   rn_c02
+   REAL(wp) ::   rn_c02r
+   REAL(wp) ::   rn_c03
+   REAL(wp) ::   rn_c04
+   REAL(wp) ::   rn_c03_sqrt2_galp
+   REAL(wp) ::   rn_sbc_mb
+   REAL(wp) ::   rn_sbc_std
+   REAL(wp) ::   rn_sbc_tke1
+   REAL(wp) ::   rn_sbc_tke2
+   REAL(wp) ::   rn_sbc_tke3
+   REAL(wp) ::   rn_sbc_psi1
+   REAL(wp) ::   rn_sbc_psi2
+   REAL(wp) ::   rn_sbc_psi3
+   REAL(wp) ::   rn_sbc_zs
+   REAL(wp) ::   zfact_tke, zfact_psi
+   REAL(wp) ::   rn_c0, rn_c2, rn_c3, rn_f6, rn_cff, rn_c_diff
+   REAL(wp) ::   rn_s0, rn_s1, rn_s2, rn_s4, rn_s5, rn_s6
+   REAL(wp) ::   rn_d0, rn_d1, rn_d2, rn_d3, rn_d4, rn_d5
+   REAL(wp) ::   rn_sc_tke, rn_sc_psi, rn_psi1, rn_psi2, rn_psi3, rn_sc_psi0
+   REAL(wp) ::   rn_psi3m, rn_psi3p, zpp, zmm, znn, epslim
+
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008) 
+   !! $Id: zdfgls.F90 1201 2008-09-24 13:24:21Z rblod $
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE zdf_gls( kt )
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE zdf_gls  ***
+      !!
+      !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity
+      !!      coefficients using a 2.5 turbulent closure scheme.
+      !!----------------------------------------------------------------------
+      USE oce,     z_elem_a  =>   ua   ! use ua as workspace
+      USE oce,     z_elem_b  =>   va   ! use va as workspace
+      USE oce,     z_elem_c  =>   ta   ! use ta as workspace
+      USE oce,     psi       =>   sa   ! use sa as workspace
+      !
+      INTEGER, INTENT(in) ::   kt ! ocean time step
+      INTEGER  ::   ji, jj, jk, ibot, ibotm1, dir  ! dummy loop arguments
+      !
+      REAL(wp) ::   zesh2, zsigpsi                 ! temporary scalars
+      REAL(wp) ::   ztx2, zty2, zup, zdown, zcof   !    -         
+      REAL(wp) ::   zratio, zrn2, zflxb, sh        !    -         -
+      REAL(wp) ::   prod, buoy, diss, zdiss, sm    !    -         -
+      REAL(wp) ::   gh, gm, shr, dif, zsqen, zav   !    -         -
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zdep   !
+      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) ::   mxlb   ! mixing length at time before
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   shear  ! vertical shear
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   eps    ! dissipation rate
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwall_psi ! Wall function for psi schmidt number in the wb case 
+                                                      ! (if ln_sigpsi.AND.ln_crban)
+      !!--------------------------------------------------------------------
+
+      IF( kt == nit000  ) CALL zdf_gls_init        ! Initialization (first time-step only)
+
+      !!--------------------------------------------------------------------
+      ! Preliminary computing
+
+      ustars2 = 0. ; ustarb2 = 0. ; psi  = 0. ; zwall_psi = 0.
+
+      ! Compute wind stress at T-points
+!CDIR NOVERRCHK
+      DO jj = 2, jpjm1
+!CDIR NOVERRCHK
+        DO ji = fs_2, fs_jpim1   ! vector opt.
+          ! 
+          ! wind stress
+          ! squared surface velocity scale
+          ztx2 = 2. * (utau(ji-1,jj  )*umask(ji-1,jj,1) + utau(ji,jj)*umask(ji,jj,1)) / &
+            &                  MAX(1., umask(ji-1,jj,1) +             umask(ji,jj,1))
+          zty2 = 2. * (vtau(ji  ,jj-1)*vmask(ji,jj-1,1) + vtau(ji,jj)*vmask(ji,jj,1)) / &
+            &                  MAX(1., vmask(ji,jj-1,1) +             vmask(ji,jj,1))
+          ustars2(ji,jj) = rn_sbc_tke2 * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)
+          !
+          ! bottom friction 
+          ztx2 = 2. * (wbotu(ji-1,jj  )*umask(ji-1,jj,1) + wbotu(ji,jj)*umask(ji,jj,1)) / &
+            &                   MAX(1., umask(ji-1,jj,1) +              umask(ji,jj,1))
+          zty2 = 2. * (wbotv(ji  ,jj-1)*vmask(ji,jj-1,1) + wbotv(ji,jj)*vmask(ji,jj,1)) / &
+            &                   MAX(1., vmask(ji,jj-1,1) +              vmask(ji,jj,1))
+          ustarb2(ji,jj) = 0.5 * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)
+        ENDDO
+      ENDDO  
+
+      ! In case of breaking surface waves mixing,
+      ! Compute surface roughness length according to Charnock formula:
+      IF (ln_crban) THEN
+        zhsro(:,:) = MAX(rn_sbc_zs * ustars2(:,:), hsro)
+      ELSE
+        zhsro(:,:) = hsro
+      ENDIF
+
+      ! 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) )   &
+                  &                            / (  fse3uw_n(ji,jj,jk)               &
+                  &                            *    fse3uw_b(ji,jj,jk) )
+               avmv(ji,jj,jk) = avmv(ji,jj,jk) * ( vn(ji,jj,jk-1) - vn(ji,jj,jk) )   &
+                  &                            * ( vb(ji,jj,jk-1) - vb(ji,jj,jk) )   &
+                  &                            / (  fse3vw_n(ji,jj,jk)               &
+                  &                            *    fse3vw_b(ji,jj,jk) )
+               eps(ji,jj,jk)  = rn_c03 * en(ji,jj,jk) * SQRT(en(ji,jj,jk)) / mxln(ji,jj,jk)
+            ENDDO
+         ENDDO
+      ENDDO
+      !
+      ! Lateral boundary conditions (avmu,avmv) (sign unchanged)
+      CALL lbc_lnk( avmu, 'U', 1. )   ;    CALL lbc_lnk( avmv, 'V', 1. )
+
+      ! Save tke at before time step
+      eb  (:,:,:) = en  (:,:,:)
+      mxlb(:,:,:) = mxln(:,:,:)
+
+      IF ( nn_clos .EQ. 0 ) THEN ! Mellor-Yamada
+         DO jk = 2, jpkm1
+            DO jj = 2, jpjm1 
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zup   = mxln(ji,jj,jk) * fsdepw(ji,jj,mbathy(ji,jj))
+                  zdown = vkarmn * fsdepw(ji,jj,jk) * (-fsdepw(ji,jj,jk)+fsdepw(ji,jj,mbathy(ji,jj)))
+                  zwall (ji,jj,jk) = ( 1. + rn_e2 * ( zup / MAX( zdown, rsmall ) )**2. ) * tmask(ji,jj,jk)
+               ENDDO
+            ENDDO
+         ENDDO
+      ENDIF
+
+      !!---------------------------------!!
+      !!   Equation to prognostic k      !!
+      !!---------------------------------!!
+      !
+      ! Now Turbulent kinetic energy (output in en)
+      ! -------------------------------
+      ! Resolution of a tridiagonal linear system by a "methode de chasse"
+      ! computation from level 2 to jpkm1  (e(1) computed after and e(jpk)=0 ).
+      ! 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
+      ! 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 + sign(0.5,shear(ji,jj,jk)+buoy)   ! dir =1(=0) if shear(ji,jj,jk)+buoy >0(<0)
+               !
+               zesh2 = dir*(shear(ji,jj,jk)+buoy)+(1-dir)*shear(ji,jj,jk)          ! production term
+               zdiss = dir*(diss/en(ji,jj,jk))   +(1-dir)*(diss-buoy)/en(ji,jj,jk) ! dissipation term
+               !
+               ! Compute a wall function from 1. to rn_sc_psi*zwall/rn_sc_psi0
+               ! Note that as long that Dirichlet boundary conditions are NOT set at the first and last levels (GOTM style)
+               ! there is no need to set a boundary condition for zwall_psi at the top and bottom boundaries.
+               ! Otherwise, this should be rn_sc_psi/rn_sc_psi0
+               IF (ln_sigpsi) THEN
+                 zsigpsi = MIN(1., zesh2/eps(ji,jj,jk)) ! 0. <= zsigpsi <= 1.
+                 zwall_psi(ji,jj,jk) = rn_sc_psi / (zsigpsi * rn_sc_psi + (1.-zsigpsi) * rn_sc_psi0 / MAX(zwall(ji,jj,jk),1.))
+               ELSE
+                 zwall_psi(ji,jj,jk) = 1.
+               ENDIF
+               !
+               ! building the matrix
+               zcof = zfact_tke * tmask(ji,jj,jk)
+               !
+               ! lower diagonal
+               z_elem_a(ji,jj,jk) = zcof * ( avm  (ji,jj,jk  ) + avm  (ji,jj,jk-1) )   &
+                  &                      / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk  ) )
+               !
+               ! upper diagonal
+               z_elem_c(ji,jj,jk) = zcof * ( avm  (ji,jj,jk+1) + avm  (ji,jj,jk  ) )   &
+                  &                      / ( fse3t(ji,jj,jk  ) * fse3w(ji,jj,jk) )
+               !
+               ! diagonal
+               z_elem_b(ji,jj,jk) = 1. - 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)
+            END DO
+         END DO
+      END DO
+      !
+      z_elem_b(:,:,jpk) = 1.
+      !
+      ! Set surface condition on zwall_psi (1 at the bottom)
+      IF (ln_sigpsi) THEN
+        DO jj = 2, jpjm1
+          DO ji = fs_2, fs_jpim1   ! vector opt.
+            zwall_psi(ji,jj,1) = rn_sc_psi/rn_sc_psi0
+          END DO
+        END DO
+      ENDIF
+
+      ! Surface boundary condition on tke
+      ! ---------------------------------
+      !
+      SELECT CASE ( nn_tkebc_surf )
+      !
+      CASE ( 0 )             ! Dirichlet case
+      !
+      IF (ln_crban) THEN     ! Wave induced mixing case
+      !                      ! en(1) = q2(1) = 0.5 * (15.8 * Ccb)^(2/3) * u*^2
+      !                      ! balance between the production and the dissipation terms including the wave effect
+      !
+      en(:,:,1) = MAX( rn_sbc_tke1 * ustars2(:,:), rn_emin )
+      z_elem_a(:,:,1) = en(:,:,1)
+      z_elem_c(:,:,1) = 0.
+      z_elem_b(:,:,1) = 1.
+      ! 
+      ! one level below
+      en(:,:,2) = MAX( rn_sbc_tke1 * ustars2(:,:) * ( (zhsro(:,:)+fsdepw(:,:,2))/zhsro(:,:) )**rn_a_sf, rn_emin )
+      z_elem_a(:,:,2) = 0.
+      z_elem_c(:,:,2) = 0.
+      z_elem_b(:,:,2) = 1.
+      !
+      ELSE                   ! No wave induced mixing case
+      !                      ! en(1) = u*^2/C0^2  &  l(1)  = K*zs
+      !                      ! balance between the production and the dissipation terms
+      !
+      en(:,:,1) = MAX( rn_c02r * ustars2(:,:), rn_emin )
+      z_elem_a(:,:,1) = en(:,:,1) 
+      z_elem_c(:,:,1) = 0.
+      z_elem_b(:,:,1) = 1.
+      !
+      ! one level below
+      en(:,:,2) = MAX( rn_c02r * ustars2(:,:), rn_emin )
+      z_elem_a(:,:,2) = 0.
+      z_elem_c(:,:,2) = 0.
+      z_elem_b(:,:,2) = 1.
+      !
+      ENDIF
+      !
+      CASE ( 1 )             ! Neumann boundary condition on d(e)/dz
+      !
+      IF (ln_crban) THEN ! Shear free case: d(e)/dz= Fw
+      !
+      ! Dirichlet conditions at k=1 (Cosmetic)
+      en(:,:,1) = MAX( rn_sbc_tke1 * ustars2(:,:), rn_emin )
+      z_elem_a(:,:,1) = en(:,:,1)
+      z_elem_c(:,:,1) = 0.
+      z_elem_b(:,:,1) = 1.
+      ! at k=2, set de/dz=Fw
+      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.           
+      zflxs(:,:) = rn_sbc_tke3 * ustars2(:,:)**1.5 * ((zhsro(:,:)+fsdept(:,:,1))/zhsro(:,:) )**(1.5*rn_a_sf)
+      en(:,:,2) = en(:,:,2) + zflxs(:,:)/fse3w(:,:,2)
+      !
+      ELSE                   ! No wave induced mixing case: d(e)/dz=0.
+      !
+      ! Dirichlet conditions at k=1 (Cosmetic)
+      en(:,:,1) = MAX( rn_c02r * ustars2(:,:), rn_emin )
+      z_elem_a(:,:,1) = en(:,:,1)
+      z_elem_c(:,:,1) = 0.
+      z_elem_b(:,:,1) = 1.
+      ! at k=2 set de/dz=0.:
+      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.
+      !
+      ENDIF
+      !
+      END SELECT
+
+      ! Bottom boundary condition on tke
+      ! --------------------------------
+      !
+      SELECT CASE ( nn_tkebc_bot )
+      !
+      CASE ( 0 )             ! Dirichlet 
+      !                      ! en(ibot) = u*^2 / Co2 and mxln(ibot) = rn_lmin
+      !                      ! Balance between the production and the dissipation terms
+      !                      
+!CDIR NOVERRCHK
+      DO jj = 2, jpjm1
+!CDIR NOVERRCHK
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            ibot   = mbathy(ji,jj)
+            ibotm1 = ibot-1
+            !
+            ! Bottom level Dirichlet condition:
+            z_elem_a(ji,jj,ibot  ) = 0.
+            z_elem_c(ji,jj,ibot  ) = 0.
+            z_elem_b(ji,jj,ibot  ) = 1.
+            en(ji,jj,ibot  ) = MAX( rn_c02r * ustarb2(ji,jj), rn_emin )
+            !
+            ! Just above last level, Dirichlet condition again
+            z_elem_a(ji,jj,ibotm1) = 0.
+            z_elem_c(ji,jj,ibotm1) = 0.
+            z_elem_b(ji,jj,ibotm1) = 1.
+            en(ji,jj,ibotm1) = MAX( rn_c02r * ustarb2(ji,jj), rn_emin ) 
+         END DO
+      END DO
+      !
+      CASE ( 1 )             ! Neumman boundary condition
+      !                      
+!CDIR NOVERRCHK
+      DO jj = 2, jpjm1
+!CDIR NOVERRCHK
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            ibot   = mbathy(ji,jj)
+            ibotm1 = ibot-1
+            !
+            ! Bottom level Dirichlet condition:
+            z_elem_a(ji,jj,ibot) = 0.
+            z_elem_c(ji,jj,ibot) = 0.
+            z_elem_b(ji,jj,ibot) = 1.
+            en(ji,jj,ibot) = MAX( rn_c02r * 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.
+         END DO
+      END DO
+      !
+      END SELECT
+
+      ! Matrix inversion (en prescribed at surface and the bottom)
+      ! ----------------------------------------------------------
+      !
+      DO jk = 2, jpkm1                             ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
+         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)
+            END DO
+         END DO
+      END DO
+      DO jk = 2, jpk                               ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
+         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)
+            END DO
+         END DO
+      END DO
+      DO jk = jpk-1, 2, -1                         ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
+         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)
+            END DO
+         END DO
+      END DO
+      !
+      ! set the minimum value of tke 
+      en(:,:,:) = MAX( en(:,:,:), rn_emin )
+      
+      !!----------------------------------------!!
+      !!   Solve prognostic equation for psi    !!
+      !!----------------------------------------!!
+
+      ! Set psi to previous time step value
+      !
+      SELECT CASE ( nn_clos )
+      !
+      CASE( 0 )               ! k-kl  (Mellor-Yamada)
+      !
+      DO jk = 2, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               psi(ji,jj,jk)  = en(ji,jj,jk) * mxln(ji,jj,jk)
+            ENDDO
+         ENDDO
+      ENDDO
+      !
+      CASE( 1 )               ! k-eps
+      !
+      DO jk = 2, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               psi(ji,jj,jk)  = eps(ji,jj,jk)
+            ENDDO
+         ENDDO
+      ENDDO
+      !
+      CASE( 2 )               ! k-w
+      !
+      DO jk = 2, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               psi(ji,jj,jk)  = SQRT(en(ji,jj,jk)) / (rn_c0 * mxln(ji,jj,jk))
+            ENDDO
+         ENDDO
+      ENDDO
+      !
+      CASE( 3 )               ! gen
+      !
+      DO jk = 2, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               psi(ji,jj,jk)  = rn_c02 * en(ji,jj,jk) * mxln(ji,jj,jk)**znn 
+            ENDDO
+         ENDDO
+      ENDDO
+      !
+      END SELECT
+      !
+      ! Now gls (output in psi)
+      ! -------------------------------
+      ! 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
+      ! 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.
+               !
+               ! psi / k
+               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 + sign(0.5,rn2(ji,jj,jk))
+               !
+               rn_psi3 = dir*rn_psi3m+(1-dir)*rn_psi3p
+               !
+               ! shear prod. - stratif. destruction
+               prod = rn_psi1 * zratio * shear(ji,jj,jk)
+               !
+               ! stratif. destruction
+               buoy = rn_psi3 * zratio * (- avt(ji,jj,jk) * rn2(ji,jj,jk))
+               !
+               ! shear prod. - stratif. destruction
+               diss = rn_psi2 * zratio * zwall(ji,jj,jk) * eps(ji,jj,jk)
+               !
+               dir = 0.5 + sign(0.5,prod+buoy)   ! dir =1(=0) if shear(ji,jj,jk)+buoy >0(<0)
+               !
+               zesh2 = dir*(prod+buoy)          +(1-dir)*prod                      ! production term
+               zdiss = dir*(diss/psi(ji,jj,jk)) +(1-dir)*(diss-buoy)/psi(ji,jj,jk) ! dissipation term
+               !                                                        
+               ! building the matrix
+               zcof = zfact_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) )   &
+                  &                      / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk  ) )
+               ! upper diagonal
+               z_elem_c(ji,jj,jk) = zcof * ( avm  (ji,jj,jk+1) + avm  (ji,jj,jk  ) )   &
+                  &                      / ( fse3t(ji,jj,jk  ) * fse3w(ji,jj,jk) )
+               ! diagonal
+               z_elem_b(ji,jj,jk) = 1. - 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)
+            END DO
+         END DO
+      END DO
+      !
+      z_elem_b(:,:,jpk) = 1.
+
+      ! Surface boundary condition on psi
+      ! ---------------------------------
+      !
+      SELECT CASE ( nn_psibc_surf )
+      !
+      CASE ( 0 )             ! Dirichlet boundary conditions
+      !
+      IF (ln_crban) THEN     ! Wave induced mixing case
+      !                      ! en(1) = q2(1) = 0.5 * (15.8 * Ccb)^(2/3) * u*^2
+      !                      ! balance between the production and the dissipation terms including the wave effect
+      !
+      zdep(:,:) = rn_l_sf * zhsro(:,:)
+      psi (:,:,1) = rn_c0**zpp * en(:,:,1)**zmm * zdep(:,:)**znn * tmask(:,:,1)
+      z_elem_a(:,:,1) = psi(:,:,1)
+      z_elem_c(:,:,1) = 0.
+      z_elem_b(:,:,1) = 1.
+      !
+      ! one level below
+      zdep(:,:) = ( (zhsro(:,:) + fsdepw(:,:,2))**(zmm*rn_a_sf+znn) )  &
+        &                       / zhsro(:,:)**(zmm*rn_a_sf)
+      psi (:,:,2) = rn_sbc_psi1 * ustars2(:,:)**zmm * zdep(:,:) * tmask(:,:,1)
+      z_elem_a(:,:,2) = 0.
+      z_elem_c(:,:,2) = 0.
+      z_elem_b(:,:,2) = 1.
+      ! 
+      ELSE                   ! No wave induced mixing case
+      !                      ! en(1) = u*^2/C0^2  &  l(1)  = K*zs
+      !                      ! balance between the production and the dissipation terms
+      !
+      zdep(:,:) = vkarmn * zhsro(:,:)
+      psi (:,:,1) = rn_c0**zpp * en(:,:,1)**zmm * zdep(:,:)**znn * tmask(:,:,1)
+      z_elem_a(:,:,1) = psi(:,:,1)
+      z_elem_c(:,:,1) = 0.
+      z_elem_b(:,:,1) = 1.
+      !
+      ! one level below
+      zdep(:,:) = vkarmn * ( zhsro(:,:) + fsdepw(:,:,2) )
+      psi (:,:,2) = rn_c0**zpp * en(:,:,1)**zmm * zdep(:,:)**znn * tmask(:,:,1)
+      z_elem_a(:,:,2) = 0.
+      z_elem_c(:,:,2) = 0.
+      z_elem_b(:,:,2) = 1.
+      !
+      ENDIF
+      !
+      CASE ( 1 )             ! Neumann boundary condition on d(psi)/dz
+      !
+      IF (ln_crban) THEN     ! Wave induced mixing case
+      !
+      zdep(:,:) = rn_l_sf * zhsro(:,:)
+      psi (:,:,1) = rn_c0**zpp * en(:,:,1)**zmm * zdep(:,:)**znn * tmask(:,:,1)
+      z_elem_a(:,:,1) = psi(:,:,1)
+      z_elem_c(:,:,1) = 0.
+      z_elem_b(:,:,1) = 1.
+      !
+      ! 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.
+      !
+      ! Set psi vertical flux at the surface:
+      zdep(:,:) = (zhsro(:,:) + fsdept(:,:,1))**(zmm*rn_a_sf+znn-1.) / zhsro(:,:)**(zmm*rn_a_sf)
+      zflxs(:,:) = rn_sbc_psi3 * ( zwall_psi(:,:,1)*avm(:,:,1) + zwall_psi(:,:,2)*avm(:,:,2) ) & 
+                 & * en(:,:,1)**zmm * zdep         
+      psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2)
+      !
+      ELSE                   ! No wave induced mixing
+      !
+      zdep(:,:) = vkarmn * zhsro(:,:)
+      psi (:,:,1) = rn_c0**zpp * en(:,:,1)**zmm * zdep(:,:)**znn * tmask(:,:,1)
+      z_elem_a(:,:,1) = psi(:,:,1)
+      z_elem_c(:,:,1) = 0.
+      z_elem_b(:,:,1) = 1.
+      !
+      ! 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(ji,jj,2) = 0.
+      !
+      ! Set psi vertical flux at the surface:
+      zdep(:,:) = zhsro(:,:) + fsdept(:,:,1)
+      zflxs(:,:) = rn_sbc_psi2 * ( avm(:,:,1) + avm(:,:,2) ) * en(:,:,1)**zmm * zdep**(znn-1.)
+      psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2)
+      !     
+      ENDIF
+      !
+      END SELECT
+
+      ! Bottom boundary condition on psi
+      ! --------------------------------
+      !
+      SELECT CASE ( nn_psibc_bot )
+      !
+      !
+      CASE ( 0 )             ! Dirichlet 
+      !                      ! en(ibot) = u*^2 / Co2 and mxln(ibot) = vkarmn * rn_hbro
+      !                      ! Balance between the production and the dissipation terms
+!CDIR NOVERRCHK
+      DO jj = 2, jpjm1
+!CDIR NOVERRCHK
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            ibot = mbathy(ji,jj)
+            ibotm1 = ibot-1
+            zdep(ji,jj) = vkarmn * rn_hbro
+            psi (ji,jj,ibot) = rn_c0**zpp * en(ji,jj,ibot)**zmm * zdep(ji,jj)**znn
+            z_elem_a(ji,jj,ibot) = 0.
+            z_elem_c(ji,jj,ibot) = 0.
+            z_elem_b(ji,jj,ibot) = 1.
+            !
+            ! Just above last level, Dirichlet condition again (GOTM like)
+            zdep(ji,jj) = vkarmn * (rn_hbro + fse3t(ji,jj,ibotm1))
+            psi (ji,jj,ibotm1) = rn_c0**zpp * en(ji,jj,ibot  )**zmm * zdep(ji,jj)**znn
+            z_elem_a(ji,jj,ibotm1) = 0.
+            z_elem_c(ji,jj,ibotm1) = 0.
+            z_elem_b(ji,jj,ibotm1) = 1.
+         END DO
+      END DO
+      !
+      !
+      CASE ( 1 )             ! Neumman boundary condition
+      !                      
+!CDIR NOVERRCHK
+      DO jj = 2, jpjm1
+!CDIR NOVERRCHK
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            ibot   = mbathy(ji,jj)
+            ibotm1 = ibot-1
+            !
+            ! Bottom level Dirichlet condition:
+            zdep(ji,jj) = vkarmn * rn_hbro
+            psi (ji,jj,ibot) = rn_c0**zpp * en(ji,jj,ibot)**zmm * zdep(ji,jj)**znn
+            !
+            z_elem_a(ji,jj,ibot) = 0.
+            z_elem_c(ji,jj,ibot) = 0.
+            z_elem_b(ji,jj,ibot) = 1.
+            !
+            ! 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.
+            !
+            ! Set psi vertical flux at the bottom:
+            zdep(ji,jj) = rn_hbro + 0.5*fse3t(ji,jj,ibotm1)
+            zflxb = rn_sbc_psi2 * ( avm(ji,jj,ibot) + avm(ji,jj,ibotm1) ) * &
+              & (0.5*(en(ji,jj,ibot)+en(ji,jj,ibotm1)))**zmm * zdep(ji,jj)**(znn-1.)
+            psi(ji,jj,ibotm1) = psi(ji,jj,ibotm1) + zflxb / fse3w(ji,jj,ibotm1)
+         END DO
+      END DO
+      !
+      END SELECT
+
+      ! Matrix inversion
+      ! ----------------
+      !
+      DO jk = 2, jpkm1                             ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
+         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)
+            END DO
+         END DO
+      END DO
+      DO jk = 2, jpk                               ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
+         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)
+            END DO
+         END DO
+      END DO
+      DO jk = jpk-1, 2, -1                         ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
+         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)
+            END DO
+         END DO
+      END DO
+
+      ! Set dissipation
+      !----------------
+
+      SELECT CASE ( nn_clos )
+      !
+      CASE( 0 )               ! k-kl  (Mellor-Yamada)
+      !
+      DO jk = 1, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+              eps(ji,jj,jk) = rn_c03 * en(ji,jj,jk) * en(ji,jj,jk) * SQRT(en(ji,jj,jk)) / psi(ji,jj,jk)
+            ENDDO
+         ENDDO
+      ENDDO
+      !
+      CASE( 1 )               ! k-eps
+      !
+      DO jk = 1, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               eps(ji,jj,jk)  = psi(ji,jj,jk)
+            ENDDO
+         ENDDO
+      ENDDO
+      !
+      CASE( 2 )               ! k-w
+      !
+      DO jk = 1, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+              eps(ji,jj,jk) = rn_c04 * en(ji,jj,jk) * psi(ji,jj,jk) 
+            ENDDO
+         ENDDO
+      ENDDO
+      !
+      CASE( 3 )               ! gen
+      !
+      DO jk = 1, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+              eps(ji,jj,jk) = rn_c0**(3.+zpp/znn) * en(ji,jj,jk)**(1.5+zmm/znn) * psi(ji,jj,jk)**(-1./znn)
+            ENDDO
+         ENDDO
+      ENDDO
+      !
+      END SELECT
+
+      ! Limit dissipation rate under stable stratification
+      ! --------------------------------------------------
+      DO jk = 1, jpkm1 ! Note that this set boundary conditions on mxln 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)  = rn_c03 * 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 )
+               mxln(ji,jj,jk) = MIN( clim_galp*SQRT( 2.*en(ji,jj,jk)/zrn2 ), mxln(ji,jj,jk) )
+            END DO
+         END DO
+      END DO 
+
+      !
+      ! Stability function and vertical viscosity and diffusivity
+      ! ---------------------------------------------------------
+      !
+      SELECT CASE ( nn_stab_func )
+      !
+      CASE ( 0 , 1 )             ! Galperin or Kantha-Clayson stability functions
+      !
+      DO jk = 2, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               ! zcof =  l²/q²
+               zcof = mxlb(ji,jj,jk)**2. / ( 2.*eb(ji,jj,jk) )
+               ! Gh = -N²l²/q²
+               gh = - rn2(ji,jj,jk) * zcof
+!               gh = MIN(gh, gh-(gh-rn_ghcri)**2./(gh + rn_gh0 - 2.0*rn_ghcri))
+               gh = MIN( gh, rn_gh0   )
+               gh = MAX( gh, rn_ghmin )
+               ! Stability functions from Kantha and Clayson (if C2=C3=0 => Galperin)
+               sh = rn_a2*( 1.-6.*rn_a1/rn_b1 ) / ( 1.-3.*rn_a2*gh*(6.*rn_a1+rn_b2*( 1.-rn_c3 ) ) )
+               sm = ( rn_b1**(-1./3.) + ( 18*rn_a1*rn_a1 + 9.*rn_a1*rn_a2*(1.-rn_c2) )*sh*gh ) / (1.-9.*rn_a1*rn_a2*gh)
+               !
+               ! Store stability function in avmu and avmv
+               avmu(ji,jj,jk) = rn_c_diff * sh * tmask(ji,jj,jk)
+               avmv(ji,jj,jk) = rn_c_diff * sm * tmask(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      !
+      CASE ( 2, 3 )               ! Canuto stability functions
+      !
+      DO jk = 2, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               ! zcof =  l²/q²
+               zcof = mxlb(ji,jj,jk)**2. / ( 2.*eb(ji,jj,jk) )
+               ! Gh = -N²l²/q²
+               gh = - rn2(ji,jj,jk) * zcof
+!               gh = MIN(gh, gh-(gh-rn_ghcri)**2./(gh + rn_gh0 - 2.0*rn_ghcri))
+               gh = MIN( gh, rn_gh0   )
+               gh = MAX( gh, rn_ghmin )
+               gh = gh*rn_f6
+               ! Gm =  M²l²/q² Shear number
+               shr = shear(ji,jj,jk)/MAX(avm(ji,jj,jk), rsmall)
+               gm = shr * zcof
+               gm = MAX(gm, 1.e-10)
+               gm = gm*rn_f6
+               gm = MIN ( (rn_d0 - rn_d1*gh + rn_d3*gh**2.)/(rn_d2-rn_d4*gh) , gm )
+               ! Stability functions from Canuto
+               rn_cff = rn_d0 - rn_d1*gh +rn_d2*gm + rn_d3*gh**2. - rn_d4*gh*gm + rn_d5*gm**2.
+               sm = (rn_s0 - rn_s1*gh + rn_s2*gm) / rn_cff
+               sh = (rn_s4 - rn_s5*gh + rn_s6*gm) / rn_cff
+               !
+               ! Store stability function in avmu and avmv
+               avmu(ji,jj,jk) = rn_c_diff * sh * tmask(ji,jj,jk)
+               avmv(ji,jj,jk) = rn_c_diff * sm * tmask(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      !
+      END SELECT
+
+      ! Boundary conditions on stability functions for momentum (Neumann):
+      ! Lines below are useless if GOTM style Dirichlet conditions are used
+      DO jj = 2, jpjm1
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            avmv(ji,jj,1) = rn_cm_sf / SQRT(2.)
+         END DO
+      END DO
+      DO jj = 2, jpjm1
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            ibot=mbathy(ji,jj)
+            ibotm1=ibot-1
+            avmv(ji,jj,ibot) = rn_c0 / SQRT(2.)
+         END DO
+      END DO
+
+      ! Compute diffusivities/viscosities
+      ! The computation below could be restrained to jk=2 to jpkm1 if GOTM style Dirichlet conditions are used
+      DO jk = 1, jpk
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zsqen           = SQRT(2.*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
+            END DO
+         END DO
+      END DO
+
+      !
+      ! Lateral boundary conditions (sign unchanged)
+      !
+      avt(:,:,1)  = 0.
+      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) = ( avm  (ji,jj,jk)*tmask(ji,jj,jk) + avm (ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk) ) * umask(ji,jj,jk)   &
+               &     / MAX( 1.,   tmask(ji,jj,jk) + tmask(ji+1,jj  ,jk) )
+               avmv(ji,jj,jk) = ( avm  (ji,jj,jk)*tmask(ji,jj,jk) + avm (ji  ,jj+1,jk)*tmask(ji  ,jj+1,jk) ) * vmask(ji,jj,jk)   &
+               &     / MAX( 1.,   tmask(ji,jj,jk) + tmask(ji  ,jj+1,jk) )
+            END DO
+         END DO
+      END DO
+
+      avmu(:,:,1) = 0.
+      avmv(:,:,1) = 0.
+
+      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 )
+      ENDIF
+      !
+   END SUBROUTINE zdf_gls
+
+   SUBROUTINE zdf_gls_init
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE zdf_gls_init  ***
+      !!                     
+      !! ** Purpose :   Initialization of the vertical eddy diffivity and 
+      !!      viscosity when 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
+      !!
+      !! ** Action  :   Increase by 1 the nstop flag is setting problem encounter
+      !!
+      !!----------------------------------------------------------------------
+      USE dynzdf_exp
+      USE trazdf_exp
+      !
+# if defined key_vectopt_memory
+      INTEGER ::   ji, jj, jk   ! dummy loop indices
+# else
+      INTEGER ::           jk   ! dummy loop indices
+# endif
+      REAL(wp)::   zcr
+      !!
+      NAMELIST/namzdf_gls/rn_emin, rn_epsmin, ln_length_lim, &
+         &            clim_galp, ln_crban, ln_sigpsi,    &
+         &            rn_crban, rn_charn,                &
+         &            nn_tkebc_surf, nn_tkebc_bot,       &
+         &            nn_psibc_surf, nn_psibc_bot,       &
+         &            nn_stab_func, nn_clos
+      !!----------------------------------------------------------
+
+      ! Read Namelist namzdf_gls
+      ! ------------------------
+      REWIND ( numnam )
+      READ   ( numnam, namzdf_gls )
+
+      ! Parameter control and print
+      ! ---------------------------
+      ! Control print
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'zdf_gls_init : gls turbulent closure scheme'
+         WRITE(numout,*) '~~~~~~~~~~~~'
+         WRITE(numout,*) '          Namelist namzdf_gls : set gls mixing parameters'
+         WRITE(numout,*) '             minimum value of en                       rn_emin  = ', rn_emin
+         WRITE(numout,*) '             minimum value of eps                    rn_epsmin  = ', rn_epsmin
+         WRITE(numout,*) '             Surface roughness (m)                         hsro = ', hsro
+         WRITE(numout,*) '             Bottom roughness (m)                        rn_hbro = ', rn_hbro
+         WRITE(numout,*) '             Limit dissipation rate under stable stratification ln_length_lim = ',ln_length_lim
+         WRITE(numout,*) '             Galperin length scale limitation coef (Standard: 0.53, Holt: 0.26) clim_galp = ', clim_galp
+         WRITE(numout,*) '             TKE Surface boundary condition       nn_tkebc_surf = ', nn_tkebc_surf
+         WRITE(numout,*) '             TKE Bottom boundary condition        nn_tkebc_bot  = ', nn_tkebc_bot
+         WRITE(numout,*) '             PSI Surface boundary condition       nn_psibc_surf = ', nn_psibc_surf
+         WRITE(numout,*) '             PSI Bottom boundary condition        nn_psibc_bot  = ', nn_psibc_bot
+         WRITE(numout,*) '             Craig and Banner scheme                  ln_crban  = ', ln_crban
+         WRITE(numout,*) '             Modify psi Schmidt number (wb case)     ln_sigpsi  = ', ln_sigpsi
+         WRITE(numout,*) '             Craig and Banner coef.                   rn_crban  = ', rn_crban
+         WRITE(numout,*) '             Charnock coef.                           rn_charn  = ', rn_charn
+         WRITE(numout,*) '             Stability functions                   nn_stab_func = ',nn_stab_func
+         WRITE(numout,*) '             Closure                                 nn_clos    = ',nn_clos
+         WRITE(numout,*)
+      ENDIF
+
+      ! Check of some namelist values
+      IF( nn_tkebc_surf < 0 .OR. nn_tkebc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_tkebc_surf is 0 or 1' )
+      IF( nn_psibc_surf < 0 .OR. nn_psibc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_psibc_surf is 0 or 1' )
+      IF( nn_tkebc_bot < 0 .OR. nn_tkebc_bot > 1 ) CALL ctl_stop( 'bad flag: nn_tkebc_bot is 0 or 1' )
+      IF( nn_psibc_bot < 0 .OR. nn_psibc_bot > 1 ) CALL ctl_stop( 'bad flag: nn_psibc_bot is 0 or 1' )
+      IF( nn_stab_func < 0 .OR. nn_stab_func > 3) CALL ctl_stop( 'bad flag: nn_stab_func is 0, 1, 2 and 3' )
+      IF( nn_clos < 0 .OR. nn_clos > 3) CALL ctl_stop( 'bad flag: nn_clos is 0, 1, 2 or 3' )
+
+      ! Initialisation of the parameters for the choosen closure
+      ! --------------------------------------------------------
+      !
+      SELECT CASE ( nn_clos )
+      !
+      CASE( 0 )               ! k-kl  (Mellor-Yamada)
+      !
+      IF(lwp) WRITE(numout,*) 'The choosen closure is k-kl closed to the classical Mellor-Yamada'
+      zpp       = 0.
+      zmm       = 1.
+      znn       = 1.
+      rn_sc_tke = 1.96
+      rn_sc_psi = 1.96
+      rn_psi1 = 0.9
+      rn_psi3p = 1.
+      rn_psi2  = 0.5
+      !
+         SELECT CASE ( nn_stab_func )
+         !
+         CASE( 0, 1 )               ! G88 or KC stability functions
+            rn_psi3m = 2.53
+         CASE( 2 )                  ! Canuto A stability functions
+            rn_psi3m = 2.38
+         CASE( 3 )                  ! Canuto B stability functions
+            rn_psi3m = 2.38         ! caution : constant not identified 
+         END SELECT
+      !
+      CASE( 1 )               ! k-eps
+      !
+      IF(lwp) WRITE(numout,*) 'The choosen closure is k-eps'
+      zpp       = 3.
+      zmm       = 1.5
+      znn       = -1.
+      rn_sc_tke = 1.
+      rn_sc_psi = 1.3  ! Schmidt number for psi
+      rn_psi1 = 1.44
+      rn_psi3p = 1.
+      rn_psi2  = 1.92
+      !
+        SELECT CASE ( nn_stab_func )
+        ! 
+        CASE( 0, 1 )               ! G88 or KC stability functions
+           rn_psi3m = -0.52
+        CASE( 2 )                  ! Canuto A stability functions
+           rn_psi3m = -0.629 
+        CASE( 3 )                  ! Canuto B stability functions
+           rn_psi3m = -0.566
+        END SELECT
+      !
+      CASE( 2 )               ! k-omega
+      !
+      IF(lwp) WRITE(numout,*) 'The choosen closure is k-omega'
+      zpp       = -1.
+      zmm       = 0.5
+      znn       = -1.
+      rn_sc_tke = 2.
+      rn_sc_psi = 2.
+      rn_psi1 = 0.555
+      rn_psi3p = 1.
+      rn_psi2  = 0.833
+      !
+        SELECT CASE ( nn_stab_func )
+        !
+        CASE( 0, 1 )               ! G88 or KC stability functions
+           rn_psi3m = -0.58
+        CASE( 2 )                  ! Canuto A stability functions
+           rn_psi3m = -0.64
+        CASE( 3 )                  ! Canuto B stability functions
+           rn_psi3m = -0.64        ! caution : constant not identified
+        END SELECT
+      !
+      CASE( 3 )               ! gen
+      !
+      IF(lwp) WRITE(numout,*) 'The choosen closure is generic'
+      zpp       = 2.
+      zmm       = 1.
+      znn       = -0.67
+      rn_sc_tke = 0.8
+      rn_sc_psi = 1.07
+      rn_psi1 = 1.
+      rn_psi3p = 1.
+      rn_psi2  = 1.22
+      !
+        SELECT CASE ( nn_stab_func )
+        !
+        CASE( 0, 1 )               ! G88 or KC stability functions
+           rn_psi3m = 0.1
+        CASE( 2 )                  ! Canuto A stability functions
+           rn_psi3m = 0.05
+        CASE( 3 )                  ! Canuto B stability functions
+           rn_psi3m = 0.05         ! caution : constant not identified
+        END SELECT
+      !
+      END SELECT
+
+      ! Initialisation of the parameters of the stability functions
+      ! -----------------------------------------------------------
+      !
+      SELECT CASE ( nn_stab_func )
+      !
+      CASE ( 0 )             ! Galperin stability functions
+      !
+      IF(lwp) WRITE(numout,*) 'Stability functions from Galperin'
+      rn_c2 = 0.
+      rn_c3 = 0.
+      rn_c_diff = 1.
+      rn_c0     = 0.5544
+      rn_cm_sf = 0.9884
+      rn_ghmin = -0.28
+      rn_gh0   = 0.0233
+      rn_ghcri = 0.02
+      !
+      CASE ( 1 )             ! Kantha-Clayson stability functions
+      !
+      IF(lwp) WRITE(numout,*) 'Stability functions from Kantha-Clayson'
+      rn_c2 = 0.7
+      rn_c3 = 0.2
+      rn_c_diff = 1.
+      rn_c0     = 0.5544
+      rn_cm_sf = 0.9884
+      rn_ghmin = -0.28
+      rn_gh0   = 0.0233
+      rn_ghcri = 0.02
+      !
+      CASE ( 2 )             ! Canuto A stability functions
+      !
+      IF(lwp) WRITE(numout,*) 'Stability functions from Canuto A'
+      rn_s0 = 1.5*rn_l1*rn_l5**2.
+      rn_s1 = -rn_l4*(rn_l6+rn_l7) + 2.*rn_l4*rn_l5*(rn_l1-(1./3.)*rn_l2-rn_l3)+1.5*rn_l1*rn_l5*rn_l8
+      rn_s2 = -(3./8.)*rn_l1*(rn_l6**2.-rn_l7**2.)
+      rn_s4 = 2.*rn_l5
+      rn_s5 = 2.*rn_l4
+      rn_s6 = (2./3.)*rn_l5*(3.*rn_l3**2.-rn_l2**2.)-0.5*rn_l5*rn_l1*(3.*rn_l3-rn_l2)+0.75*rn_l1*(rn_l6-rn_l7)
+      rn_d0 = 3*rn_l5**2.
+      rn_d1 = rn_l5*(7.*rn_l4+3.*rn_l8)
+      rn_d2 = rn_l5**2.*(3.*rn_l3**2.-rn_l2**2.)-0.75*(rn_l6**2.-rn_l7**2.)
+      rn_d3 = rn_l4*(4.*rn_l4+3.*rn_l8)
+      rn_d4 = rn_l4*(rn_l2*rn_l6-3.*rn_l3*rn_l7-rn_l5*(rn_l2**2.-rn_l3**2.))+rn_l5*rn_l8*(3.*rn_l3**2.-rn_l2**2.)
+      rn_d5 = 0.25*(rn_l2**2.-3.*rn_l3**2.)*(rn_l6**2.-rn_l7**2.)
+      rn_c0 = 0.5268
+      rn_f6 = 8. / (rn_c0**6.)
+      rn_c_diff = SQRT(2.)/(rn_c0**3.)
+      rn_cm_sf = 0.7310
+      rn_ghmin = -0.28
+      rn_gh0   = 0.0329
+      rn_ghcri = 0.03
+      !
+      CASE ( 3 )             ! Canuto B stability functions
+      !
+      IF(lwp) WRITE(numout,*) 'Stability functions from Canuto B'
+      rn_s0 = 1.5*rn_m1*rn_m5**2.
+      rn_s1 = -rn_m4*(rn_m6+rn_m7) + 2.*rn_m4*rn_m5*(rn_m1-(1./3.)*rn_m2-rn_m3)+1.5*rn_m1*rn_m5*rn_m8
+      rn_s2 = -(3./8.)*rn_m1*(rn_m6**2.-rn_m7**2.)
+      rn_s4 = 2.*rn_m5
+      rn_s5 = 2.*rn_m4
+      rn_s6 = (2./3.)*rn_m5*(3.*rn_m3**2.-rn_m2**2.)-0.5*rn_m5*rn_m1*(3.*rn_m3-rn_m2)+0.75*rn_m1*(rn_m6-rn_m7)
+      rn_d0 = 3*rn_m5**2.
+      rn_d1 = rn_m5*(7.*rn_m4+3.*rn_m8)
+      rn_d2 = rn_m5**2.*(3.*rn_m3**2.-rn_m2**2.)-0.75*(rn_m6**2.-rn_m7**2.)
+      rn_d3 = rn_m4*(4.*rn_m4+3.*rn_m8)
+      rn_d4 = rn_m4*(rn_m2*rn_m6-3.*rn_m3*rn_m7-rn_m5*(rn_m2**2.-rn_m3**2.))+rn_m5*rn_m8*(3.*rn_m3**2.-rn_m2**2.)
+      rn_d5 = 0.25*(rn_m2**2.-3.*rn_m3**2.)*(rn_m6**2.-rn_m7**2.)
+      rn_c0 = 0.5268  !!       rn_c0 = 0.5540 (Warner ...) to verify !
+      rn_f6 = 8. / (rn_c0**6.)
+      rn_c_diff = SQRT(2.)/(rn_c0**3.)
+      rn_cm_sf = 0.7470
+      rn_ghmin = -0.28
+      rn_gh0   = 0.0444
+      rn_ghcri = 0.0414
+      !
+      END SELECT
+    
+      ! Set Schmidt number for psi diffusion
+      ! In the wave breaking case
+      ! See equation 13 of Carniel et al, Ocean modelling, 30, 225-239, 2009
+      ! or equation (17) of Burchard, JPO, 31, 3133-3145, 2001
+      IF ((ln_sigpsi).AND.(ln_crban)) THEN
+         zcr = SQRT(1.5*rn_sc_tke) * rn_cm_sf /vkarmn
+         rn_sc_psi0 = vkarmn**2/(rn_psi2*rn_cm_sf**2) *              & 
+        &             (  znn**2 - 4./3.*zcr*znn*zmm - 1./3.*zcr*znn  &
+        &              + 2./9.*zmm*zcr**2 + 4./9.*zcr**2*zmm**2)                                 
+      ELSE
+         rn_sc_psi0 = rn_sc_psi
+      ENDIF
+ 
+      ! Shear free turbulence parameters:
+      !
+      rn_a_sf  = -4.*znn*SQRT(rn_sc_tke) / ( (1.+4.*zmm)*SQRT(rn_sc_tke) &
+               &                              - SQRT(rn_sc_tke + 24.*rn_sc_psi0*rn_psi2 ) )
+      rn_l_sf  = rn_c0 * SQRT(rn_c0/rn_cm_sf) * SQRT( ( (1. + 4.*zmm + 8.*zmm**2)*rn_sc_tke       &
+               &                                          + 12. * rn_sc_psi0*rn_psi2 - (1. + 4.*zmm) &
+               &                                          *SQRT(rn_sc_tke*(rn_sc_tke                &
+               &                                                + 24.*rn_sc_psi0*rn_psi2)) )         &
+               &                                          /(12.*znn**2.) &
+               &                                       )
+
+      ! Control print
+      !
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'Limit values'
+         WRITE(numout,*) '~~~~~~~~~~~~'
+         WRITE(numout,*) 'Parameter  m = ',zmm
+         WRITE(numout,*) 'Parameter  n = ',znn
+         WRITE(numout,*) 'Parameter  p = ',zpp
+         WRITE(numout,*) 'rn_psi1   = ',rn_psi1
+         WRITE(numout,*) 'rn_psi2   = ',rn_psi2
+         WRITE(numout,*) 'rn_psi3m  = ',rn_psi3m
+         WRITE(numout,*) 'rn_psi3p  = ',rn_psi3p
+         WRITE(numout,*) 'rn_sc_tke = ',rn_sc_tke
+         WRITE(numout,*) 'rn_sc_psi = ',rn_sc_psi
+         WRITE(numout,*) 'rn_sc_psi0 = ',rn_sc_psi0
+         WRITE(numout,*) 'rn_c0     = ',rn_c0
+         WRITE(numout,*)
+         WRITE(numout,*) 'Shear free turbulence parameters:'
+         WRITE(numout,*) 'rn_cm_sf  = ',rn_cm_sf
+         WRITE(numout,*) 'rn_a_sf   = ',rn_a_sf
+         WRITE(numout,*) 'rn_l_sf   = ',rn_l_sf
+         WRITE(numout,*)
+      ENDIF
+
+      ! Constants initialization
+      rn_c02r = 1. / rn_c0**2.
+      rn_c02  = rn_c0**2._wp
+      rn_c03  = rn_c0**3._wp
+      rn_c04  = rn_c0**4._wp
+      rn_c03_sqrt2_galp = rn_c03 / SQRT(2._wp) / clim_galp
+      rn_sbc_mb   = 0.5 * (15.8*rn_crban)**(2./3.) ! Surf. bound. cond. from Mellor and Blumberg
+      rn_sbc_std  = 3.75                           ! Surf. bound. cond. standard (prod=diss)
+      rn_sbc_tke1 = (-rn_sc_tke*rn_crban/(rn_cm_sf*rn_a_sf*rn_l_sf))**(2./3.) ! k_eps = 53.  Dirichlet + Wave breaking 
+      rn_sbc_tke2 = 0.5 / rau0
+      rn_sbc_tke3 = rdt * rn_crban                                                         ! Neumann + Wave breaking
+      rn_sbc_zs   = rn_charn/grav                                                          ! Charnock formula
+      rn_sbc_psi1 = rn_c0**zpp * rn_sbc_tke1**zmm * rn_l_sf**znn                           ! Dirichlet + Wave breaking
+      rn_sbc_psi2 = -0.5 * rdt * rn_c0**zpp * znn * vkarmn**znn / rn_sc_psi                   ! Neumann + NO Wave breaking 
+      rn_sbc_psi3 = -0.5 * rdt * rn_c0**zpp * rn_l_sf**znn / rn_sc_psi  * (znn + zmm*rn_a_sf) ! Neumann + Wave breaking
+      zfact_tke = -0.5 / rn_sc_tke * rdt               ! Cst used for the Diffusion term of tke
+      zfact_psi = -0.5 / rn_sc_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
+      !                               !* read or initialize all required files 
+      CALL gls_rst( nit000, 'READ' )
+      !
+   END SUBROUTINE zdf_gls_init
+
+   SUBROUTINE gls_rst( kt, cdrw )
+     !!---------------------------------------------------------------------
+     !!                   ***  ROUTINE ts_rst  ***
+     !!                     
+     !! ** Purpose :   Read or write TKE file (en) in restart file
+     !!
+     !! ** Method  :   use of IOM library
+     !!                if the restart does not contain TKE, en is either 
+     !!                set to rn_emin or recomputed (nn_igls/=0)
+     !!----------------------------------------------------------------------
+     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, id7, id8
+     INTEGER ::   ji, jj, ikbu, ikbv, ikbum1, ikbvm1
+     REAL(wp)::   cbx, cby
+     !!----------------------------------------------------------------------
+     !
+     IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise 
+        !                                   ! ---------------
+        IF( ln_rstart ) THEN                   !* Read the restart file
+           id1 = iom_varid( numror, 'en'   , ldstop = .FALSE. )
+           id2 = iom_varid( numror, 'avt'  , ldstop = .FALSE. )
+           id3 = iom_varid( numror, 'avm'  , ldstop = .FALSE. )
+           id4 = iom_varid( numror, 'avmu' , ldstop = .FALSE. )
+           id5 = iom_varid( numror, 'avmv' , ldstop = .FALSE. )
+           id6 = iom_varid( numror, 'mxln' , ldstop = .FALSE. )
+           id7 = iom_varid( numror, 'wbotu', ldstop = .FALSE. )
+           id8 = iom_varid( numror, 'wbotv', ldstop = .FALSE. )
+           !
+           IF( MIN( id1, id2, id3, id4, id5, id6, id7, id8 ) > 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, 'wbotu' , wbotu  )
+              CALL iom_get( numror, jpdom_autoglo, 'wbotv' , wbotv  )
+           ELSE                        
+              IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without gls scheme, en and mxln computed by iterative loop'
+              IF(lwp) WRITE(numout,*) ' ===>>>> : The bottom stresses are estimated' 
+              en  (:,:,:) = rn_emin
+              mxln(:,:,:) = 0.001        
+              ! Initialize bottom stresses
+              DO jj = 2, jpjm1
+                DO ji = fs_2, fs_jpim1   ! vector opt.
+                  ikbu   = MIN( mbathy(ji+1,jj  ), mbathy(ji,jj) )
+                  ikbum1 = MAX( ikbu-1, 1 )
+                  ikbv   = MIN( mbathy(ji,jj+1), mbathy(ji,jj) )
+                  ikbvm1 = MAX( ikbv-1, 1 )
+                  cbx = avmu(ji,jj,ikbu) / fse3uw(ji,jj,ikbu)
+                  cby = avmv(ji,jj,ikbv) / fse3vw(ji,jj,ikbv)
+                  wbotu(ji,jj) = -cbx * un(ji,jj,ikbum1)*umask(ji,jj,1)
+                  wbotv(ji,jj) = -cby * vn(ji,jj,ikbvm1)*vmask(ji,jj,1)
+                END DO
+              END DO
+              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,*) ' ===>>>> : The bottom stresses are estimated'
+           en  (:,:,:) = rn_emin
+           mxln(:,:,:) = 0.001       
+           ! Initialize bottom stresses
+           DO jj = 2, jpjm1
+             DO ji = fs_2, fs_jpim1   ! vector opt.
+               ikbu   = MIN( mbathy(ji+1,jj  ), mbathy(ji,jj) )
+               ikbum1 = MAX( ikbu-1, 1 )
+               ikbv   = MIN( mbathy(ji,jj+1), mbathy(ji,jj) )
+               ikbvm1 = MAX( ikbv-1, 1 )
+               cbx = avmu(ji,jj,ikbu) / fse3uw(ji,jj,ikbu)
+               cby = avmv(ji,jj,ikbv) / fse3vw(ji,jj,ikbv)
+               wbotu(ji,jj) = -cbx * un(ji,jj,ikbum1)*umask(ji,jj,1)
+               wbotv(ji,jj) = -cby * vn(ji,jj,ikbvm1)*vmask(ji,jj,1)
+             END DO
+           END DO
+        ENDIF
+        !
+     ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
+        !                                   ! -------------------
+        IF(lwp) WRITE(numout,*) '---- gls-rst ----'
+        CALL iom_rstput( kt, nitrst, numrow, 'en'   , en    )
+        CALL iom_rstput( kt, nitrst, numrow, 'avt'  , avt   )
+        CALL iom_rstput( kt, nitrst, numrow, 'avm'  , avm   )
+        CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu  )
+        CALL iom_rstput( kt, nitrst, numrow, 'avmv' , avmv  )
+        CALL iom_rstput( kt, nitrst, numrow, 'mxln' , mxln  )
+        !
+     ENDIF
+     !
+   END SUBROUTINE gls_rst
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Dummy module :                                        NO TKE scheme
+   !!----------------------------------------------------------------------
+   LOGICAL, PUBLIC, PARAMETER ::   lk_zdfgls = .FALSE.   !: TKE flag
+CONTAINS
+   SUBROUTINE zdf_gls( kt )          ! Empty routine
+      WRITE(*,*) 'zdf_gls: You should not have seen this print! error?', kt
+   END SUBROUTINE zdf_gls
+#endif
+
+   !!======================================================================
+END MODULE zdfgls
Index: branches/devmercator2010_1/NEMO/OPA_SRC/ZDF/zdfini.F90
===================================================================
--- branches/devmercator2010_1/NEMO/OPA_SRC/ZDF/zdfini.F90	(revision 2131)
+++ branches/devmercator2010_1/NEMO/OPA_SRC/ZDF/zdfini.F90	(revision 2132)
@@ -20,4 +20,5 @@
    USE zdftke_old      ! TKE vertical mixing  (old scheme)
    USE zdftke          ! TKE vertical mixing
+   USE zdfgls          ! GLS vertical mixing
    USE zdfkpp          ! KPP vertical mixing          
    USE zdfddm          ! double diffusion mixing      
@@ -106,4 +107,8 @@
          ioptio = ioptio+1
       ENDIF
+      IF( lk_zdfgls ) THEN
+         IF(lwp) WRITE(numout,*) '      GLS dependent eddy coefficients'
+         ioptio = ioptio+1
+      ENDIF
       IF( lk_zdfkpp ) THEN
          IF(lwp) WRITE(numout,*) '      KPP dependent eddy coefficients'
@@ -128,4 +133,7 @@
          IF(lwp) WRITE(numout,*) '      use the 1.5 turbulent closure'
       ENDIF
+      IF( lk_zdfgls ) THEN
+         IF(lwp) WRITE(numout,*) '      use the GLS closure scheme'
+      ENDIF
       IF( lk_zdfkpp ) THEN
          IF(lwp) WRITE(numout,*) '      use the KPP closure scheme'
@@ -136,6 +144,6 @@
       ENDIF
       IF ( ioptio > 1 .AND. .NOT. lk_esopa )   CALL ctl_stop( ' chose between ln_zdfnpc and ln_zdfevd' )
-      IF( ioptio == 0 .AND. .NOT.( lk_zdftke_old .OR. lk_zdftke .OR. lk_zdfkpp ) ) &
-         CALL ctl_stop( ' except for TKE sor KPP physics, a convection scheme is', &
+      IF( ioptio == 0 .AND. .NOT.( lk_zdftke_old .OR. lk_zdftke .OR. lk_zdfgls .OR. lk_zdfkpp ) ) &
+         CALL ctl_stop( ' except for TKE, GLS or KPP physics, a convection scheme is', &
          &              ' required: ln_zdfevd or ln_zdfnpc logicals' )
 
Index: branches/devmercator2010_1/NEMO/OPA_SRC/step.F90
===================================================================
--- branches/devmercator2010_1/NEMO/OPA_SRC/step.F90	(revision 2131)
+++ branches/devmercator2010_1/NEMO/OPA_SRC/step.F90	(revision 2132)
@@ -90,4 +90,5 @@
    USE zdftke_old      ! old TKE vertical mixing      (zdf_tke_old routine)
    USE zdftke          ! TKE vertical mixing              (zdf_tke routine)
+   USE zdfgls          ! GLS vertical mixing              (zdf_gls routine)  
    USE zdfkpp          ! KPP vertical mixing              (zdf_kpp routine)
    USE zdfddm          ! double diffusion mixing          (zdf_ddm routine)
@@ -204,4 +205,5 @@
       IF( lk_zdftke_old )   CALL zdf_tke_old( kstp )  ! TKE closure scheme for Kz (old scheme)
       IF( lk_zdftke     )   CALL zdf_tke    ( kstp )  ! TKE closure scheme for Kz
+      IF( lk_zdfgls     )   CALL zdf_gls    ( kstp )  ! GLS closure scheme for Kz
       IF( lk_zdfkpp     )   CALL zdf_kpp    ( kstp )  ! KPP closure scheme for Kz
       IF( lk_zdfcst     )   THEN                      ! Constant Kz (reset avt, avm[uv] to the background value)
@@ -223,4 +225,6 @@
                                                       ! write tke information in the restart file
       IF( lrst_oce .AND. lk_zdftke )   CALL tke_rst( kstp, 'WRITE' )
+                                                      ! write gls information in the restart file
+      IF( lrst_oce .AND. lk_zdfgls )   CALL gls_rst( kstp, 'WRITE' )
       !
       !  LATERAL  PHYSICS 
