Index: /branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90
===================================================================
--- /branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90	(revision 2396)
+++ /branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90	(revision 2397)
@@ -5,6 +5,6 @@
    !!                 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
    !!----------------------------------------------------------------------
 #if defined key_zdfgls   ||   defined key_esopa
@@ -34,5 +34,5 @@
 
    PUBLIC   zdf_gls        ! routine called in step module
-   PUBLIC   zdf_gls_init   ! routine called in step module
+   PUBLIC   zdf_gls_init   ! routine called in opa module
    PUBLIC   gls_rst        ! routine called in step module
 
@@ -46,8 +46,6 @@
    !                                         !!! ** Namelist  namzdf_gls  **
    LOGICAL  ::   ln_crban      = .FALSE.      ! =T use Craig and Banner scheme
-   LOGICAL  ::   ln_length_lim = .FALSE.      ! use limit on the dissipation rate under stable stratif. (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)
+   LOGICAL  ::   ln_length_lim = .FALSE.      ! use limit on the dissipation rate under stable stratification (Galperin et al. 1988)
+   LOGICAL  ::   ln_sigpsi     = .FALSE.      ! Activate Burchard (2003) modification for k-eps closure & wave breaking mixing
    INTEGER  ::   nn_tkebc_surf = 0            ! TKE surface boundary condition (=0/1)
    INTEGER  ::   nn_tkebc_bot  = 0            ! TKE bottom boundary condition (=0/1)
@@ -57,14 +55,16 @@
    INTEGER  ::   nn_clos       = 0            ! closure 0/1/2/3 MY82/k-eps/k-w/gen
    REAL(wp) ::   rn_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_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)
    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) ::   rcm_sf        = 0.73_wp        ! Shear free turbulence parameters
-   REAL(wp) ::   ra_sf         = -2.0_wp        ! Must be negative -2 < ra_sf < -1 
-   REAL(wp) ::   rl_sf         =  0.2_wp        ! 0 <rl_sf<vkarmn    
-   REAL(wp) ::   rghmin        = -0.28
-   REAL(wp) ::   rgh0          = 0.0329
-   REAL(wp) ::   rghcri        = 0.03
+   REAL(wp) ::   hsro          =  0.003_wp    ! Minimum surface roughness
+   REAL(wp) ::   rcm_sf        =  0.73_wp     ! Shear free turbulence parameters
+   REAL(wp) ::   ra_sf         = -2.0_wp      ! Must be negative -2 < ra_sf < -1 
+   REAL(wp) ::   rl_sf         =  0.2_wp      ! 0 <rl_sf<vkarmn    
+   REAL(wp) ::   rghmin        = -0.28_wp
+   REAL(wp) ::   rgh0          =  0.0329_wp
+   REAL(wp) ::   rghcri        =  0.03_wp
    REAL(wp) ::   ra1           =  0.92_wp
    REAL(wp) ::   ra2           =  0.74_wp
@@ -88,24 +88,15 @@
    REAL(wp) ::   rm7           =  0.0_wp
    REAL(wp) ::   rm8           =  0.318_wp
-   REAL(wp) ::   rc02
-   REAL(wp) ::   rc02r
-   REAL(wp) ::   rc03
-   REAL(wp) ::   rc04
-   REAL(wp) ::   rc03_sqrt2_galp
-   REAL(wp) ::   rsbc_mb
-   REAL(wp) ::   rsbc_std
-   REAL(wp) ::   rsbc_tke1
-   REAL(wp) ::   rsbc_tke2
-   REAL(wp) ::   rsbc_tke3
-   REAL(wp) ::   rsbc_psi1
-   REAL(wp) ::   rsbc_psi2
-   REAL(wp) ::   rsbc_psi3
-   REAL(wp) ::   rsbc_zs
-   REAL(wp) ::   rfact_tke, rfact_psi
-   REAL(wp) ::   rc0, rc2, rc3, rf6, rcff, rc_diff
-   REAL(wp) ::   rs0, rs1, rs2, rs4, rs5, rs6
-   REAL(wp) ::   rd0, rd1, rd2, rd3, rd4, rd5
-   REAL(wp) ::   rsc_tke, rsc_psi, rpsi1, rpsi2, rpsi3, rsc_psi0
-   REAL(wp) ::   rpsi3m, rpsi3p, rpp, rmm, rnn
+   
+   REAL(wp) ::   rc02, rc02r, rc03, rc04                          ! coefficients deduced from above parameters
+   REAL(wp) ::   rc03_sqrt2_galp                                  !     -           -           -        -
+   REAL(wp) ::   rsbc_tke1, rsbc_tke2, rsbc_tke3, rfact_tke       !     -           -           -        -
+   REAL(wp) ::   rsbc_psi1, rsbc_psi2, rsbc_psi3, rfact_psi       !     -           -           -        -
+   REAL(wp) ::   rsbc_mb  , rsbc_std , rsbc_zs                    !     -           -           -        -
+   REAL(wp) ::   rc0, rc2, rc3, rf6, rcff, rc_diff                !     -           -           -        -
+   REAL(wp) ::   rs0, rs1, rs2, rs4, rs5, rs6                     !     -           -           -        -
+   REAL(wp) ::   rd0, rd1, rd2, rd3, rd4, rd5                     !     -           -           -        -
+   REAL(wp) ::   rsc_tke, rsc_psi, rpsi1, rpsi2, rpsi3, rsc_psi0  !     -           -           -        -
+   REAL(wp) ::   rpsi3m, rpsi3p, rpp, rmm, rnn                    !     -           -           -        -
 
    !! * Substitutions
@@ -124,5 +115,5 @@
       !!
       !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity
-      !!      coefficients using a 2.5 turbulent closure scheme.
+      !!              coefficients using the GLS turbulent closure scheme.
       !!----------------------------------------------------------------------
       USE oce,     z_elem_a  =>   ua   ! use ua as workspace
@@ -133,24 +124,22 @@
       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)
+      REAL(wp) ::   zesh2, zsigpsi, zcoef, zex1, zex2   ! local 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 use in the wb case (ln_sigpsi.AND.ln_crban=T)
       !!--------------------------------------------------------------------
 
       ! Preliminary computing
 
-      ustars2 = 0. ; ustarb2 = 0. ; psi  = 0. ; zwall_psi = 0.
+      ustars2 = 0._wp   ;   ustarb2 = 0._wp   ;   psi  = 0._wp   ;   zwall_psi = 0._wp
 
       ! Compute wind stress at T-points
@@ -162,25 +151,19 @@
           ! 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) = rsbc_tke2 * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)
+          ustars2(ji,jj) = rau0r * taum(ji,jj) * 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  
+          ztx2 = ( 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 = ( 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) = SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)
+        END DO
+      END DO  
 
       ! In case of breaking surface waves mixing,
       ! Compute surface roughness length according to Charnock formula:
-      IF (ln_crban) THEN
-        zhsro(:,:) = MAX(rsbc_zs * ustars2(:,:), hsro)
-      ELSE
-        zhsro(:,:) = hsro
+      IF( ln_crban ) THEN   ;   zhsro(:,:) = MAX(rsbc_zs * ustars2(:,:), hsro)
+      ELSE                  ;   zhsro(:,:) = hsro
       ENDIF
 
@@ -198,10 +181,10 @@
                   &                            *    fse3vw_b(ji,jj,jk) )
                eps(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT(en(ji,jj,jk)) / mxln(ji,jj,jk)
-            ENDDO
-         ENDDO
-      ENDDO
+            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. )
+      CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )
 
       ! Save tke at before time step
@@ -209,14 +192,15 @@
       mxlb(:,:,:) = mxln(:,:,:)
 
-      IF ( nn_clos .EQ. 0 ) THEN ! Mellor-Yamada
+      IF( nn_clos == 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. + re2 * ( zup / MAX( zdown, rsmall ) )**2. ) * tmask(ji,jj,jk)
-               ENDDO
-            ENDDO
-         ENDDO
+                  zdown = vkarmn * fsdepw(ji,jj,jk) * ( -fsdepw(ji,jj,jk) + fsdepw(ji,jj,mbathy(ji,jj)) )
+                  zcoef = ( zup / MAX( zdown, rsmall ) )
+                  zwall (ji,jj,jk) = ( 1._wp + re2 * zcoef*zcoef ) * tmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
       ENDIF
 
@@ -248,8 +232,8 @@
                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
+               dir = 0.5_wp + SIGN( 0.5_wp, shear(ji,jj,jk) + buoy )   ! dir =1(=0) if shear(ji,jj,jk)+buoy >0(<0)
+               !
+               zesh2 = dir*(shear(ji,jj,jk)+buoy)+(1._wp-dir)*shear(ji,jj,jk)          ! production term
+               zdiss = dir*(diss/en(ji,jj,jk))   +(1._wp-dir)*(diss-buoy)/en(ji,jj,jk) ! dissipation term
                !
                ! Compute a wall function from 1. to rsc_psi*zwall/rsc_psi0
@@ -257,9 +241,9 @@
                ! there is no need to set a boundary condition for zwall_psi at the top and bottom boundaries.
                ! Otherwise, this should be rsc_psi/rsc_psi0
-               IF (ln_sigpsi) THEN
-                 zsigpsi = MIN(1., zesh2/eps(ji,jj,jk)) ! 0. <= zsigpsi <= 1.
-                 zwall_psi(ji,jj,jk) = rsc_psi / (zsigpsi * rsc_psi + (1.-zsigpsi) * rsc_psi0 / MAX(zwall(ji,jj,jk),1.))
+               IF( ln_sigpsi ) THEN
+                  zsigpsi = MIN( 1._wp, zesh2 / eps(ji,jj,jk) )     ! 0. <= zsigpsi <= 1.
+                  zwall_psi(ji,jj,jk) = rsc_psi / (  zsigpsi * rsc_psi + (1._wp-zsigpsi) * rsc_psi0 / MAX( zwall(ji,jj,jk), 1._wp )  )
                ELSE
-                 zwall_psi(ji,jj,jk) = 1.
+                  zwall_psi(ji,jj,jk) = 1._wp
                ENDIF
                !
@@ -276,6 +260,6 @@
                !
                ! 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) 
+               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
@@ -285,13 +269,14 @@
       END DO
       !
-      z_elem_b(:,:,jpk) = 1.
+      z_elem_b(:,:,jpk) = 1._wp
       !
       ! 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) = rsc_psi/rsc_psi0
-          END DO
-        END DO
+      IF( ln_sigpsi ) THEN
+         zcoef = rsc_psi / rsc_psi0
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zwall_psi(ji,jj,1) = zcoef
+            END DO
+         END DO
       ENDIF
 
@@ -302,65 +287,63 @@
       !
       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( rsbc_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( rsbc_tke1 * ustars2(:,:) * ( (zhsro(:,:)+fsdepw(:,:,2))/zhsro(:,:) )**ra_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( rc02r * 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( rc02r * ustars2(:,:), rn_emin )
-      z_elem_a(:,:,2) = 0.
-      z_elem_c(:,:,2) = 0.
-      z_elem_b(:,:,2) = 1.
-      !
-      ENDIF
-      !
+         !
+         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( rsbc_tke1 * ustars2(:,:), rn_emin )
+            z_elem_a(:,:,1) = en(:,:,1)
+            z_elem_c(:,:,1) = 0._wp
+            z_elem_b(:,:,1) = 1._wp
+            ! 
+            ! one level below
+            en(:,:,2) = MAX( rsbc_tke1 * ustars2(:,:) * ( (zhsro(:,:)+fsdepw(:,:,2))/zhsro(:,:) )**ra_sf, rn_emin )
+            z_elem_a(:,:,2) = 0._wp
+            z_elem_c(:,:,2) = 0._wp
+            z_elem_b(:,:,2) = 1._wp
+            !
+         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( rc02r * ustars2(:,:), rn_emin )
+            z_elem_a(:,:,1) = en(:,:,1) 
+            z_elem_c(:,:,1) = 0._wp
+            z_elem_b(:,:,1) = 1._wp
+            !
+            ! one level below
+            en(:,:,2) = MAX( rc02r * ustars2(:,:), rn_emin )
+            z_elem_a(:,:,2) = 0._wp
+            z_elem_c(:,:,2) = 0._wp
+            z_elem_b(:,:,2) = 1._wp
+            !
+         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( rsbc_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(:,:) = rsbc_tke3 * ustars2(:,:)**1.5 * ((zhsro(:,:)+fsdept(:,:,1))/zhsro(:,:) )**(1.5*ra_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( rc02r * 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
-      !
+         !
+         IF (ln_crban) THEN ! Shear free case: d(e)/dz= Fw
+            !
+            ! Dirichlet conditions at k=1 (Cosmetic)
+            en(:,:,1) = MAX( rsbc_tke1 * ustars2(:,:), rn_emin )
+            z_elem_a(:,:,1) = en(:,:,1)
+            z_elem_c(:,:,1) = 0._wp
+            z_elem_b(:,:,1) = 1._wp
+            ! 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._wp        
+            zflxs(:,:) = rsbc_tke3 * ustars2(:,:)**1.5_wp * ( (zhsro(:,:)+fsdept(:,:,1) ) / zhsro(:,:) )**(1.5*ra_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( rc02r * ustars2(:,:), rn_emin )
+            z_elem_a(:,:,1) = en(:,:,1)
+            z_elem_c(:,:,1) = 0._wp
+            z_elem_b(:,:,1) = 1._wp
+            ! 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._wp
+            !
+         ENDIF
+         !
       END SELECT
 
@@ -371,49 +354,48 @@
       !
       CASE ( 0 )             ! Dirichlet 
-      !                      ! en(ibot) = u*^2 / Co2 and mxln(ibot) = rn_lmin
-      !                      ! Balance between the production and the dissipation terms
-      !                      
+         !                      ! en(ibot) = u*^2 / Co2 and mxln(ibot) = rn_lmin
+         !                      ! Balance between the production and the dissipation terms
 !CDIR NOVERRCHK
-      DO jj = 2, jpjm1
+         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( rc02r * 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( rc02r * ustarb2(ji,jj), rn_emin ) 
-         END DO
-      END DO
-      !
+            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._wp
+               z_elem_c(ji,jj,ibot  ) = 0._wp
+               z_elem_b(ji,jj,ibot  ) = 1._wp
+               en(ji,jj,ibot  ) = MAX( rc02r * ustarb2(ji,jj), rn_emin )
+               !
+               ! Just above last level, Dirichlet condition again
+               z_elem_a(ji,jj,ibotm1) = 0._wp
+               z_elem_c(ji,jj,ibotm1) = 0._wp
+               z_elem_b(ji,jj,ibotm1) = 1._wp
+               en(ji,jj,ibotm1) = MAX( rc02r * ustarb2(ji,jj), rn_emin ) 
+            END DO
+         END DO
+         !
       CASE ( 1 )             ! Neumman boundary condition
-      !                      
+         !                      
 !CDIR NOVERRCHK
-      DO jj = 2, jpjm1
+         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( rc02r * ustarb2(ji,jj), rn_emin )
-            !
-            ! Just above last level: Neumann condition
-            z_elem_b(ji,jj,ibotm1) = z_elem_b(ji,jj,ibotm1) + z_elem_c(ji,jj,ibotm1) ! Remove z_elem_c from z_elem_b
-            z_elem_c(ji,jj,ibotm1) = 0.
-         END DO
-      END DO
-      !
+            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._wp
+               z_elem_c(ji,jj,ibot) = 0._wp
+               z_elem_b(ji,jj,ibot) = 1._wp
+               en(ji,jj,ibot) = MAX( rc02r * ustarb2(ji,jj), rn_emin )
+               !
+               ! Just above last level: Neumann condition
+               z_elem_b(ji,jj,ibotm1) = z_elem_b(ji,jj,ibotm1) + z_elem_c(ji,jj,ibotm1)   ! Remove z_elem_c from z_elem_b
+               z_elem_c(ji,jj,ibotm1) = 0._wp
+            END DO
+         END DO
+         !
       END SELECT
 
@@ -442,6 +424,5 @@
          END DO
       END DO
-      !
-      ! set the minimum value of tke 
+      !                                            ! set the minimum value of tke 
       en(:,:,:) = MAX( en(:,:,:), rn_emin )
       
@@ -455,43 +436,39 @@
       !
       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
-      !
+         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)
+               END DO
+            END DO
+         END DO
+         !
       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
-      !
+         DO jk = 2, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  psi(ji,jj,jk)  = eps(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+         !
       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)) / (rc0 * 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)  = rc02 * en(ji,jj,jk) * mxln(ji,jj,jk)**rnn 
-            ENDDO
-         ENDDO
-      ENDDO
-      !
+         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) ) / ( rc0 * mxln(ji,jj,jk) )
+               END DO
+            END DO
+         END DO
+         !
+      CASE( 3 )               ! generic
+         DO jk = 2, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  psi(ji,jj,jk)  = rc02 * en(ji,jj,jk) * mxln(ji,jj,jk)**rnn 
+               END DO
+            END DO
+         END DO
+         !
       END SELECT
       !
@@ -511,7 +488,7 @@
                !
                ! 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))
-               !
-               rpsi3 = dir*rpsi3m+(1-dir)*rpsi3p
+               dir = 0.5_wp + SIGN( 0.5_wp, rn2(ji,jj,jk) )
+               !
+               rpsi3 = dir * rpsi3m + ( 1._wp - dir ) * rpsi3p
                !
                ! shear prod. - stratif. destruction
@@ -519,13 +496,13 @@
                !
                ! stratif. destruction
-               buoy = rpsi3 * zratio * (- avt(ji,jj,jk) * rn2(ji,jj,jk))
+               buoy = rpsi3 * zratio * (- avt(ji,jj,jk) * rn2(ji,jj,jk) )
                !
                ! shear prod. - stratif. destruction
                diss = rpsi2 * zratio * zwall(ji,jj,jk) * eps(ji,jj,jk)
                !
-               dir = 0.5 + 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
+               dir = 0.5_wp + SIGN( 0.5_wp, prod + buoy )   ! dir =1(=0) if shear(ji,jj,jk)+buoy >0(<0)
+               !
+               zesh2 = dir * ( prod + buoy )          + (1._wp - dir ) * prod                        ! production term
+               zdiss = dir * ( diss / psi(ji,jj,jk) ) + (1._wp - dir ) * (diss-buoy) / psi(ji,jj,jk) ! dissipation term
                !                                                        
                ! building the matrix
@@ -538,6 +515,6 @@
                   &                      / ( 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)
+               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
@@ -547,5 +524,5 @@
       END DO
       !
-      z_elem_b(:,:,jpk) = 1.
+      z_elem_b(:,:,jpk) = 1._wp
 
       ! Surface boundary condition on psi
@@ -555,81 +532,81 @@
       !
       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(:,:) = rl_sf * zhsro(:,:)
-      psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * 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))**(rmm*ra_sf+rnn) )  &
-        &                       / zhsro(:,:)**(rmm*ra_sf)
-      psi (:,:,2) = rsbc_psi1 * ustars2(:,:)**rmm * 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) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * 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) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
-      z_elem_a(:,:,2) = 0.
-      z_elem_c(:,:,2) = 0.
-      z_elem_b(:,:,2) = 1.
-      !
-      ENDIF
-      !
+         !
+         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(:,:) = rl_sf * zhsro(:,:)
+            psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
+            z_elem_a(:,:,1) = psi(:,:,1)
+            z_elem_c(:,:,1) = 0._wp
+            z_elem_b(:,:,1) = 1._wp
+            !
+            ! one level below
+            zex1 = (rmm*ra_sf+rnn)
+            zex2 = (rmm*ra_sf)
+            zdep(:,:) = ( (zhsro(:,:) + fsdepw(:,:,2))**zex1 ) / zhsro(:,:)**zex2
+            psi (:,:,2) = rsbc_psi1 * ustars2(:,:)**rmm * zdep(:,:) * tmask(:,:,1)
+            z_elem_a(:,:,2) = 0._wp
+            z_elem_c(:,:,2) = 0._wp
+            z_elem_b(:,:,2) = 1._wp
+            ! 
+         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) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
+            z_elem_a(:,:,1) = psi(:,:,1)
+            z_elem_c(:,:,1) = 0._wp
+            z_elem_b(:,:,1) = 1._wp
+            !
+            ! one level below
+            zdep(:,:) = vkarmn * ( zhsro(:,:) + fsdepw(:,:,2) )
+            psi (:,:,2) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
+            z_elem_a(:,:,2) = 0._wp
+            z_elem_c(:,:,2) = 0._wp
+            z_elem_b(:,:,2) = 1.
+            !
+         ENDIF
+         !
       CASE ( 1 )             ! Neumann boundary condition on d(psi)/dz
-      !
-      IF (ln_crban) THEN     ! Wave induced mixing case
-      !
-      zdep(:,:) = rl_sf * zhsro(:,:)
-      psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * 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))**(rmm*ra_sf+rnn-1.) / zhsro(:,:)**(rmm*ra_sf)
-      zflxs(:,:) = rsbc_psi3 * ( zwall_psi(:,:,1)*avm(:,:,1) + zwall_psi(:,:,2)*avm(:,:,2) ) & 
-                 & * en(:,:,1)**rmm * zdep         
-      psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2)
-      !
+         !
+         IF( ln_crban ) THEN     ! Wave induced mixing case
+            !
+            zdep(:,:) = rl_sf * zhsro(:,:)
+            psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
+            z_elem_a(:,:,1) = psi(:,:,1)
+            z_elem_c(:,:,1) = 0._wp
+            z_elem_b(:,:,1) = 1._wp
+            !
+            ! Neumann condition at k=2
+            z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b
+            z_elem_a(:,:,2) = 0._wp
+            !
+            ! Set psi vertical flux at the surface:
+            zdep(:,:) = (zhsro(:,:) + fsdept(:,:,1))**(rmm*ra_sf+rnn-1._wp) / zhsro(:,:)**(rmm*ra_sf)
+            zflxs(:,:) = rsbc_psi3 * ( zwall_psi(:,:,1)*avm(:,:,1) + zwall_psi(:,:,2)*avm(:,:,2) ) & 
+               &                   * en(:,:,1)**rmm * zdep         
+            psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2)
+            !
       ELSE                   ! No wave induced mixing
-      !
-      zdep(:,:) = vkarmn * zhsro(:,:)
-      psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * 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(:,:) = rsbc_psi2 * ( avm(:,:,1) + avm(:,:,2) ) * en(:,:,1)**rmm * zdep**(rnn-1.)
-      psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2)
-      !     
-      ENDIF
-      !
+            !
+            zdep(:,:) = vkarmn * zhsro(:,:)
+            psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
+            z_elem_a(:,:,1) = psi(:,:,1)
+            z_elem_c(:,:,1) = 0._wp
+            z_elem_b(:,:,1) = 1._wp
+            !
+            ! Neumann condition at k=2
+            z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b
+            z_elem_a(ji,jj,2) = 0._wp
+            !
+            ! Set psi vertical flux at the surface:
+            zdep(:,:)  = zhsro(:,:) + fsdept(:,:,1)
+            zflxs(:,:) = rsbc_psi2 * ( avm(:,:,1) + avm(:,:,2) ) * en(:,:,1)**rmm * zdep**(rnn-1._wp)
+            psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2)
+            !     
+         ENDIF
+         !
       END SELECT
 
@@ -639,59 +616,57 @@
       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
+         !                      ! 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
+         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) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn
-            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) = rc0**rpp * en(ji,jj,ibot  )**rmm * zdep(ji,jj)**rnn
-            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
-      !
-      !
+            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) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn
+               z_elem_a(ji,jj,ibot) = 0._wp
+               z_elem_c(ji,jj,ibot) = 0._wp
+               z_elem_b(ji,jj,ibot) = 1._wp
+               !
+               ! Just above last level, Dirichlet condition again (GOTM like)
+               zdep(ji,jj) = vkarmn * (rn_hbro + fse3t(ji,jj,ibotm1))
+               psi (ji,jj,ibotm1) = rc0**rpp * en(ji,jj,ibot  )**rmm * zdep(ji,jj)**rnn
+               z_elem_a(ji,jj,ibotm1) = 0._wp
+               z_elem_c(ji,jj,ibotm1) = 0._wp
+               z_elem_b(ji,jj,ibotm1) = 1._wp
+            END DO
+         END DO
+         !
       CASE ( 1 )             ! Neumman boundary condition
-      !                      
+         !                      
 !CDIR NOVERRCHK
-      DO jj = 2, jpjm1
+         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) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn
-            !
-            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 = rsbc_psi2 * ( avm(ji,jj,ibot) + avm(ji,jj,ibotm1) ) * &
-              & (0.5*(en(ji,jj,ibot)+en(ji,jj,ibotm1)))**rmm * zdep(ji,jj)**(rnn-1.)
-            psi(ji,jj,ibotm1) = psi(ji,jj,ibotm1) + zflxb / fse3w(ji,jj,ibotm1)
-         END DO
-      END DO
-      !
+            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) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn
+               !
+               z_elem_a(ji,jj,ibot) = 0._wp
+               z_elem_c(ji,jj,ibot) = 0._wp
+               z_elem_b(ji,jj,ibot) = 1._wp
+               !
+               ! 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_wp*fse3t(ji,jj,ibotm1)
+               zflxb = rsbc_psi2 * ( avm(ji,jj,ibot) + avm(ji,jj,ibotm1) )   &
+                  &  * (0.5_wp*(en(ji,jj,ibot)+en(ji,jj,ibotm1)))**rmm * zdep(ji,jj)**(rnn-1._wp)
+               psi(ji,jj,ibotm1) = psi(ji,jj,ibotm1) + zflxb / fse3w(ji,jj,ibotm1)
+            END DO
+         END DO
+         !
       END SELECT
 
@@ -727,43 +702,42 @@
       !
       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) = rc03 * en(ji,jj,jk) * en(ji,jj,jk) * SQRT(en(ji,jj,jk)) / psi(ji,jj,jk)
-            ENDDO
-         ENDDO
-      ENDDO
-      !
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / psi(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+         !
       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
-      !
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  eps(ji,jj,jk) = psi(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+         !
       CASE( 2 )               ! k-w
-      !
-      DO jk = 1, jpkm1
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-              eps(ji,jj,jk) = rc04 * 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) = rc0**(3.+rpp/rnn) * en(ji,jj,jk)**(1.5+rmm/rnn) * psi(ji,jj,jk)**(-1./rnn)
-            ENDDO
-         ENDDO
-      ENDDO
-      !
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  eps(ji,jj,jk) = rc04 * en(ji,jj,jk) * psi(ji,jj,jk) 
+               END DO
+            END DO
+         END DO
+         !
+      CASE( 3 )               ! generic
+         zcoef = rc0**( 3._wp  + rpp/rnn )
+         zex1  =      ( 1.5_wp + rmm/rnn )
+         zex2  = -1._wp / rnn
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  eps(ji,jj,jk) = zcoef * en(ji,jj,jk)**zex1 * psi(ji,jj,jk)**zex2
+               END DO
+            END DO
+         END DO
+         !
       END SELECT
 
@@ -775,8 +749,8 @@
                ! 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)
+               mxln(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 )
-               mxln(ji,jj,jk) = MIN( rn_clim_galp*SQRT( 2.*en(ji,jj,jk)/zrn2 ), mxln(ji,jj,jk) )
+               mxln(ji,jj,jk) = MIN(  rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), mxln(ji,jj,jk)  )
             END DO
          END DO
@@ -790,69 +764,67 @@
       !
       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, rgh0   )
-               gh = MAX( gh, rghmin )
-               ! Stability functions from Kantha and Clayson (if C2=C3=0 => Galperin)
-               sh = ra2*( 1.-6.*ra1/rb1 ) / ( 1.-3.*ra2*gh*(6.*ra1+rb2*( 1.-rc3 ) ) )
-               sm = ( rb1**(-1./3.) + ( 18*ra1*ra1 + 9.*ra1*ra2*(1.-rc2) )*sh*gh ) / (1.-9.*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)
-            END DO
-         END DO
-      END DO
-      !
+         DO jk = 2, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  ! zcof =  l²/q²
+                  zcof = mxlb(ji,jj,jk) * mxlb(ji,jj,jk) / ( 2._wp*eb(ji,jj,jk) )
+                  ! Gh = -N²l²/q²
+                  gh = - rn2(ji,jj,jk) * zcof
+                  gh = MIN( gh, rgh0   )
+                  gh = MAX( gh, rghmin )
+                  ! Stability functions from Kantha and Clayson (if C2=C3=0 => Galperin)
+                  sh = ra2*( 1._wp-6._wp*ra1/rb1 ) / ( 1.-3.*ra2*gh*(6.*ra1+rb2*( 1._wp-rc3 ) ) )
+                  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)
+               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, rgh0   )
-               gh = MAX( gh, rghmin )
-               gh = gh*rf6
-               ! 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*rf6
-               gm = MIN ( (rd0 - rd1*gh + rd3*gh**2.)/(rd2-rd4*gh) , gm )
-               ! Stability functions from Canuto
-               rcff = rd0 - rd1*gh +rd2*gm + rd3*gh**2. - rd4*gh*gm + rd5*gm**2.
-               sm = (rs0 - rs1*gh + rs2*gm) / rcff
-               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)
-            END DO
-         END DO
-      END DO
-      !
+         DO jk = 2, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  ! zcof =  l²/q²
+                  zcof = mxlb(ji,jj,jk)*mxlb(ji,jj,jk) / ( 2._wp * eb(ji,jj,jk) )
+                  ! Gh = -N²l²/q²
+                  gh = - rn2(ji,jj,jk) * zcof
+                  gh = MIN( gh, rgh0   )
+                  gh = MAX( gh, rghmin )
+                  gh = gh * rf6
+                  ! Gm =  M²l²/q² Shear number
+                  shr = shear(ji,jj,jk) / MAX( avm(ji,jj,jk), rsmall )
+                  gm = MAX( shr * zcof , 1.e-10 )
+                  gm = gm * rf6
+                  gm = MIN ( (rd0 - rd1*gh + rd3*gh*gh) / (rd2-rd4*gh) , gm )
+                  ! Stability functions from Canuto
+                  rcff = rd0 - rd1*gh +rd2*gm + rd3*gh*gh - rd4*gh*gm + rd5*gm*gm
+                  sm = (rs0 - rs1*gh + rs2*gm) / rcff
+                  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)
+               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
+      zcoef = rcm_sf / SQRT( 2._wp )
       DO jj = 2, jpjm1
          DO ji = fs_2, fs_jpim1   ! vector opt.
-            avmv(ji,jj,1) = rcm_sf / SQRT(2.)
-         END DO
-      END DO
+            avmv(ji,jj,1) = zcoef
+         END DO
+      END DO
+      zcoef = rc0 / SQRT( 2._wp )
       DO jj = 2, jpjm1
          DO ji = fs_2, fs_jpim1   ! vector opt.
-            ibot=mbathy(ji,jj)
-            ibotm1=ibot-1
-            avmv(ji,jj,ibot) = rc0 / SQRT(2.)
+            ibot   = mbathy(ji,jj)
+            avmv(ji,jj,ibot) = zcoef
          END DO
       END DO
@@ -863,17 +835,15 @@
          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
-
+               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
+            END DO
+         END DO
+      END DO
       !
       ! Lateral boundary conditions (sign unchanged)
-      !
-      avt(:,:,1)  = 0.
+      avt(:,:,1)  = 0._wp
       CALL lbc_lnk( avm, 'W', 1. )   ;   CALL lbc_lnk( avt, 'W', 1. )
 
@@ -881,16 +851,11 @@
          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
+               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
@@ -918,4 +883,7 @@
       !!
       !!----------------------------------------------------------------------
+      USE dynzdf_exp
+      USE trazdf_exp
+      !
       INTEGER ::   jk    ! dummy loop indices
       REAL(wp)::   zcr   ! local scalar
@@ -929,43 +897,38 @@
       !!----------------------------------------------------------
 
-      ! Read Namelist namzdf_gls
-      ! ------------------------
-      REWIND ( numnam )
+      REWIND ( numnam )                !* Read Namelist namzdf_gls
       READ   ( numnam, namzdf_gls )
 
-      ! Parameter control and print
-      ! ---------------------------
-      ! Control print
-      IF(lwp) THEN
+      IF(lwp) THEN                     !* Control print
          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) rn_clim_galp = ', rn_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,*) '   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 stratif.  ln_length_lim = ', ln_length_lim
+         WRITE(numout,*) '      Galperin limit (Standard: 0.53, Holt: 0.26)   rn_clim_galp  = ', rn_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 coefficient                  rn_crban       = ', rn_crban
+         WRITE(numout,*) '      Charnock coefficient                          rn_charn       = ', rn_charn
+         WRITE(numout,*) '      Stability functions                           nn_stab_func   = ', nn_stab_func
+         WRITE(numout,*) '      Type of closure                               nn_clos        = ', nn_clos
          WRITE(numout,*)
       ENDIF
 
-      ! Check of some namelist values
+      !                                !* 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' )
+      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
@@ -975,91 +938,75 @@
       !
       CASE( 0 )               ! k-kl  (Mellor-Yamada)
-      !
-      IF(lwp) WRITE(numout,*) 'The choosen closure is k-kl closed to the classical Mellor-Yamada'
-      rpp       = 0.
-      rmm       = 1.
-      rnn       = 1.
-      rsc_tke = 1.96
-      rsc_psi = 1.96
-      rpsi1 = 0.9
-      rpsi3p = 1.
-      rpsi2  = 0.5
-      !
+         !
+         IF(lwp) WRITE(numout,*) 'The choosen closure is k-kl closed to the classical Mellor-Yamada'
+         rpp     = 0._wp
+         rmm     = 1._wp
+         rnn     = 1._wp
+         rsc_tke = 1.96_wp
+         rsc_psi = 1.96_wp
+         rpsi1   = 0.9_wp
+         rpsi3p  = 1._wp
+         rpsi2   = 0.5_wp
+         !
          SELECT CASE ( nn_stab_func )
-         !
-         CASE( 0, 1 )               ! G88 or KC stability functions
-            rpsi3m = 2.53
-         CASE( 2 )                  ! Canuto A stability functions
-            rpsi3m = 2.38
-         CASE( 3 )                  ! Canuto B stability functions
-            rpsi3m = 2.38         ! caution : constant not identified 
+         CASE( 0, 1 )   ;   rpsi3m = 2.53_wp       ! G88 or KC stability functions
+         CASE( 2 )      ;   rpsi3m = 2.38_wp       ! Canuto A stability functions
+         CASE( 3 )      ;   rpsi3m = 2.38          ! Canuto B stability functions (caution : constant not identified)
          END SELECT
-      !
+         !
       CASE( 1 )               ! k-eps
-      !
-      IF(lwp) WRITE(numout,*) 'The choosen closure is k-eps'
-      rpp       = 3.
-      rmm       = 1.5
-      rnn       = -1.
-      rsc_tke = 1.
-      rsc_psi = 1.3  ! Schmidt number for psi
-      rpsi1 = 1.44
-      rpsi3p = 1.
-      rpsi2  = 1.92
-      !
-        SELECT CASE ( nn_stab_func )
-        ! 
-        CASE( 0, 1 )               ! G88 or KC stability functions
-           rpsi3m = -0.52
-        CASE( 2 )                  ! Canuto A stability functions
-           rpsi3m = -0.629 
-        CASE( 3 )                  ! Canuto B stability functions
-           rpsi3m = -0.566
-        END SELECT
-      !
+         !
+         IF(lwp) WRITE(numout,*) 'The choosen closure is k-eps'
+         rpp     =  3._wp
+         rmm     =  1.5_wp
+         rnn     = -1._wp
+         rsc_tke =  1._wp
+         rsc_psi =  1.3_wp  ! Schmidt number for psi
+         rpsi1   =  1.44_wp
+         rpsi3p  =  1._wp
+         rpsi2   =  1.92_wp
+         !
+         SELECT CASE ( nn_stab_func )
+         CASE( 0, 1 )   ;   rpsi3m = -0.52_wp      ! G88 or KC stability functions
+         CASE( 2 )      ;   rpsi3m = -0.629_wp     ! Canuto A stability functions
+         CASE( 3 )      ;   rpsi3m = -0.566        ! Canuto B stability functions
+         END SELECT
+         !
       CASE( 2 )               ! k-omega
-      !
-      IF(lwp) WRITE(numout,*) 'The choosen closure is k-omega'
-      rpp       = -1.
-      rmm       = 0.5
-      rnn       = -1.
-      rsc_tke = 2.
-      rsc_psi = 2.
-      rpsi1 = 0.555
-      rpsi3p = 1.
-      rpsi2  = 0.833
-      !
-        SELECT CASE ( nn_stab_func )
-        !
-        CASE( 0, 1 )               ! G88 or KC stability functions
-           rpsi3m = -0.58
-        CASE( 2 )                  ! Canuto A stability functions
-           rpsi3m = -0.64
-        CASE( 3 )                  ! Canuto B stability functions
-           rpsi3m = -0.64        ! caution : constant not identified
-        END SELECT
-      !
-      CASE( 3 )               ! gen
-      !
-      IF(lwp) WRITE(numout,*) 'The choosen closure is generic'
-      rpp       = 2.
-      rmm       = 1.
-      rnn       = -0.67
-      rsc_tke = 0.8
-      rsc_psi = 1.07
-      rpsi1 = 1.
-      rpsi3p = 1.
-      rpsi2  = 1.22
-      !
-        SELECT CASE ( nn_stab_func )
-        !
-        CASE( 0, 1 )               ! G88 or KC stability functions
-           rpsi3m = 0.1
-        CASE( 2 )                  ! Canuto A stability functions
-           rpsi3m = 0.05
-        CASE( 3 )                  ! Canuto B stability functions
-           rpsi3m = 0.05         ! caution : constant not identified
-        END SELECT
-      !
+         !
+         IF(lwp) WRITE(numout,*) 'The choosen closure is k-omega'
+         rpp     = -1._wp
+         rmm     =  0.5_wp
+         rnn     = -1._wp
+         rsc_tke =  2._wp
+         rsc_psi =  2._wp
+         rpsi1   =  0.555_wp
+         rpsi3p  =  1._wp
+         rpsi2   =  0.833_wp
+         !
+         SELECT CASE ( nn_stab_func )
+         CASE( 0, 1 )   ;   rpsi3m = -0.58_wp       ! G88 or KC stability functions
+         CASE( 2 )      ;   rpsi3m = -0.64_wp       ! Canuto A stability functions
+         CASE( 3 )      ;   rpsi3m = -0.64_wp       ! Canuto B stability functions caution : constant not identified)
+         END SELECT
+         !
+      CASE( 3 )               ! generic
+         !
+         IF(lwp) WRITE(numout,*) 'The choosen closure is generic'
+         rpp     = 2._wp
+         rmm     = 1._wp
+         rnn     = -0.67_wp
+         rsc_tke = 0.8_wp
+         rsc_psi = 1.07_wp
+         rpsi1   = 1._wp
+         rpsi3p  = 1._wp
+         rpsi2   = 1.22_wp
+         !
+         SELECT CASE ( nn_stab_func )
+         CASE( 0, 1 )   ;   rpsi3m = 0.1_wp         ! G88 or KC stability functions
+         CASE( 2 )      ;   rpsi3m = 0.05_wp        ! Canuto A stability functions
+         CASE( 3 )      ;   rpsi3m = 0.05_wp        ! Canuto B stability functions caution : constant not identified)
+         END SELECT
+         !
       END SELECT
 
@@ -1070,84 +1017,84 @@
       !
       CASE ( 0 )             ! Galperin stability functions
-      !
-      IF(lwp) WRITE(numout,*) 'Stability functions from Galperin'
-      rc2 = 0.
-      rc3 = 0.
-      rc_diff = 1.
-      rc0     = 0.5544
-      rcm_sf = 0.9884
-      rghmin = -0.28
-      rgh0   = 0.0233
-      rghcri = 0.02
-      !
+         !
+         IF(lwp) WRITE(numout,*) 'Stability functions from Galperin'
+         rc2     =  0._wp
+         rc3     =  0._wp
+         rc_diff =  1._wp
+         rc0     =  0.5544_wp
+         rcm_sf  =  0.9884_wp
+         rghmin  = -0.28_wp
+         rgh0    =  0.0233_wp
+         rghcri  =  0.02_wp
+         !
       CASE ( 1 )             ! Kantha-Clayson stability functions
-      !
-      IF(lwp) WRITE(numout,*) 'Stability functions from Kantha-Clayson'
-      rc2 = 0.7
-      rc3 = 0.2
-      rc_diff = 1.
-      rc0     = 0.5544
-      rcm_sf = 0.9884
-      rghmin = -0.28
-      rgh0   = 0.0233
-      rghcri = 0.02
-      !
+         !
+         IF(lwp) WRITE(numout,*) 'Stability functions from Kantha-Clayson'
+         rc2     =  0.7_wp
+         rc3     =  0.2_wp
+         rc_diff =  1._wp
+         rc0     =  0.5544_wp
+         rcm_sf  =  0.9884_wp
+         rghmin  = -0.28_wp
+         rgh0    =  0.0233_wp
+         rghcri  =  0.02_wp
+         !
       CASE ( 2 )             ! Canuto A stability functions
-      !
-      IF(lwp) WRITE(numout,*) 'Stability functions from Canuto A'
-      rs0 = 1.5*rl1*rl5**2.
-      rs1 = -rl4*(rl6+rl7) + 2.*rl4*rl5*(rl1-(1./3.)*rl2-rl3)+1.5*rl1*rl5*rl8
-      rs2 = -(3./8.)*rl1*(rl6**2.-rl7**2.)
-      rs4 = 2.*rl5
-      rs5 = 2.*rl4
-      rs6 = (2./3.)*rl5*(3.*rl3**2.-rl2**2.)-0.5*rl5*rl1*(3.*rl3-rl2)+0.75*rl1*(rl6-rl7)
-      rd0 = 3*rl5**2.
-      rd1 = rl5*(7.*rl4+3.*rl8)
-      rd2 = rl5**2.*(3.*rl3**2.-rl2**2.)-0.75*(rl6**2.-rl7**2.)
-      rd3 = rl4*(4.*rl4+3.*rl8)
-      rd4 = rl4*(rl2*rl6-3.*rl3*rl7-rl5*(rl2**2.-rl3**2.))+rl5*rl8*(3.*rl3**2.-rl2**2.)
-      rd5 = 0.25*(rl2**2.-3.*rl3**2.)*(rl6**2.-rl7**2.)
-      rc0 = 0.5268
-      rf6 = 8. / (rc0**6.)
-      rc_diff = SQRT(2.)/(rc0**3.)
-      rcm_sf = 0.7310
-      rghmin = -0.28
-      rgh0   = 0.0329
-      rghcri = 0.03
-      !
+         !
+         IF(lwp) WRITE(numout,*) 'Stability functions from Canuto A'
+         rs0 = 1.5_wp * rl1 * rl5*rl5
+         rs1 = -rl4*(rl6+rl7) + 2._wp*rl4*rl5*(rl1-(1._wp/3._wp)*rl2-rl3) + 1.5_wp*rl1*rl5*rl8
+         rs2 = -(3._wp/8._wp) * rl1*(rl6*rl6-rl7*rl7)
+         rs4 = 2._wp * rl5
+         rs5 = 2._wp * rl4
+         rs6 = (2._wp/3._wp) * rl5 * ( 3._wp*rl3*rl3 - rl2*rl2 ) - 0.5_wp * rl5*rl1 * (3._wp*rl3-rl2)   &
+            &                                                    + 0.75_wp * rl1 * ( rl6 - rl7 )
+         rd0 = 3._wp * rl5*rl5
+         rd1 = rl5 * ( 7._wp*rl4 + 3._wp*rl8 )
+         rd2 = rl5*rl5 * ( 3._wp*rl3*rl3 - rl2*rl2 ) - 0.75_wp*(rl6*rl6 - rl7*rl7 )
+         rd3 = rl4 * ( 4._wp*rl4 + 3._wp*rl8)
+         rd4 = rl4 * ( rl2 * rl6 - 3._wp*rl3*rl7 - rl5*(rl2*rl2 - rl3*rl3 ) ) + rl5*rl8 * ( 3._wp*rl3*rl3 - rl2*rl2 )
+         rd5 = 0.25_wp * ( rl2*rl2 - 3._wp *rl3*rl3 ) * ( rl6*rl6 - rl7*rl7 )
+         rc0 = 0.5268_wp
+         rf6 = 8._wp / (rc0**6._wp)
+         rc_diff = SQRT(2._wp) / (rc0**3._wp)
+         rcm_sf  =  0.7310_wp
+         rghmin  = -0.28_wp
+         rgh0    =  0.0329_wp
+         rghcri  =  0.03_wp
+         !
       CASE ( 3 )             ! Canuto B stability functions
-      !
-      IF(lwp) WRITE(numout,*) 'Stability functions from Canuto B'
-      rs0 = 1.5*rm1*rm5**2.
-      rs1 = -rm4*(rm6+rm7) + 2.*rm4*rm5*(rm1-(1./3.)*rm2-rm3)+1.5*rm1*rm5*rm8
-      rs2 = -(3./8.)*rm1*(rm6**2.-rm7**2.)
-      rs4 = 2.*rm5
-      rs5 = 2.*rm4
-      rs6 = (2./3.)*rm5*(3.*rm3**2.-rm2**2.)-0.5*rm5*rm1*(3.*rm3-rm2)+0.75*rm1*(rm6-rm7)
-      rd0 = 3*rm5**2.
-      rd1 = rm5*(7.*rm4+3.*rm8)
-      rd2 = rm5**2.*(3.*rm3**2.-rm2**2.)-0.75*(rm6**2.-rm7**2.)
-      rd3 = rm4*(4.*rm4+3.*rm8)
-      rd4 = rm4*(rm2*rm6-3.*rm3*rm7-rm5*(rm2**2.-rm3**2.))+rm5*rm8*(3.*rm3**2.-rm2**2.)
-      rd5 = 0.25*(rm2**2.-3.*rm3**2.)*(rm6**2.-rm7**2.)
-      rc0 = 0.5268  !!       rc0 = 0.5540 (Warner ...) to verify !
-      rf6 = 8. / (rc0**6.)
-      rc_diff = SQRT(2.)/(rc0**3.)
-      rcm_sf = 0.7470
-      rghmin = -0.28
-      rgh0   = 0.0444
-      rghcri = 0.0414
-      !
+         !
+         IF(lwp) WRITE(numout,*) 'Stability functions from Canuto B'
+         rs0 = 1.5_wp * rm1 * rm5*rm5
+         rs1 = -rm4 * (rm6+rm7) + 2._wp * rm4*rm5*(rm1-(1._wp/3._wp)*rm2-rm3) + 1.5_wp * rm1*rm5*rm8
+         rs2 = -(3._wp/8._wp) * rm1 * (rm6*rm6-rm7*rm7 )
+         rs4 = 2._wp * rm5
+         rs5 = 2._wp * rm4
+         rs6 = (2._wp/3._wp) * rm5 * (3._wp*rm3*rm3-rm2*rm2) - 0.5_wp * rm5*rm1*(3._wp*rm3-rm2) + 0.75_wp * rm1*(rm6-rm7)
+         rd0 = 3._wp * rm5*rm5
+         rd1 = rm5 * (7._wp*rm4 + 3._wp*rm8)
+         rd2 = rm5*rm5 * (3._wp*rm3*rm3 - rm2*rm2) - 0.75_wp * (rm6*rm6 - rm7*rm7)
+         rd3 = rm4 * ( 4._wp*rm4 + 3._wp*rm8 )
+         rd4 = rm4 * ( rm2*rm6 -3._wp*rm3*rm7 - rm5*(rm2*rm2 - rm3*rm3) ) + rm5 * rm8 * ( 3._wp*rm3*rm3 - rm2*rm2 )
+         rd5 = 0.25_wp * ( rm2*rm2 - 3._wp*rm3*rm3 ) * ( rm6*rm6 - rm7*rm7 )
+         rc0 = 0.5268_wp            !!       rc0 = 0.5540_wp (Warner ...) to verify !
+         rf6 = 8._wp / ( rc0**6._wp )
+         rc_diff = SQRT(2._wp)/(rc0**3.)
+         rcm_sf  =  0.7470_wp
+         rghmin  = -0.28_wp
+         rgh0    =  0.0444_wp
+         rghcri  =  0.0414_wp
+         !
       END SELECT
     
-      ! Set Schmidt number for psi diffusion
-      ! In the wave breaking case
+      ! 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*rsc_tke) * rcm_sf /vkarmn
-         rsc_psi0 = vkarmn**2/(rpsi2*rcm_sf**2) *              & 
-        &             (  rnn**2 - 4./3.*zcr*rnn*rmm - 1./3.*zcr*rnn  &
-        &              + 2./9.*rmm*zcr**2 + 4./9.*zcr**2*rmm**2)                                 
+      IF( ln_sigpsi .AND. ln_crban ) THEN
+         zcr = SQRT( 1.5_wp*rsc_tke ) * rcm_sf / vkarmn
+         rsc_psi0 = vkarmn*vkarmn / ( rpsi2 * rcm_sf*rcm_sf )                       & 
+        &         * ( rnn*rnn - 4._wp/3._wp * zcr*rnn*rmm - 1._wp/3._wp * zcr*rnn   &
+        &           + 2._wp/9._wp * rmm * zcr*zcr + 4._wp/9._wp * zcr*zcr * rmm*rmm )                                 
       ELSE
          rsc_psi0 = rsc_psi
@@ -1156,16 +1103,14 @@
       ! Shear free turbulence parameters:
       !
-      ra_sf  = -4.*rnn*SQRT(rsc_tke) / ( (1.+4.*rmm)*SQRT(rsc_tke) &
-               &                              - SQRT(rsc_tke + 24.*rsc_psi0*rpsi2 ) )
-      rl_sf  = rc0 * SQRT(rc0/rcm_sf) * SQRT( ( (1. + 4.*rmm + 8.*rmm**2)*rsc_tke       &
-               &                                          + 12. * rsc_psi0*rpsi2 - (1. + 4.*rmm) &
-               &                                          *SQRT(rsc_tke*(rsc_tke                &
-               &                                                + 24.*rsc_psi0*rpsi2)) )         &
-               &                                          /(12.*rnn**2.) &
-               &                                       )
-
-      ! Control print
-      !
-      IF(lwp) THEN
+      ra_sf  = -4._wp * rnn * SQRT( rsc_tke ) / ( (1._wp+4._wp*rmm) * SQRT( rsc_tke )   &
+         &                                      - SQRT(rsc_tke + 24._wp*rsc_psi0*rpsi2 ) )
+      rl_sf  = rc0 * SQRT( rc0 / rcm_sf )                                                                   &
+         &         * SQRT(  (  (1._wp + 4._wp*rmm + 8._wp*rmm*rmm) * rsc_tke                                &
+         &                   + 12._wp * rsc_psi0 * rpsi2                                                    &
+         &                   - (1._wp + 4._wp*rmm) * SQRT( rsc_tke*(rsc_tke+ 24._wp*rsc_psi0*rpsi2) )  )    &
+         &                / ( 12._wp*rnn*rnn )                                                              )
+
+      !
+      IF(lwp) THEN      ! Control print
          WRITE(numout,*)
          WRITE(numout,*) 'Limit values'
@@ -1191,25 +1136,24 @@
 
       ! Constants initialization
-      rc02r = 1. / rc0**2.
-      rc02  = rc0**2._wp
-      rc03  = rc0**3._wp
-      rc04  = rc0**4._wp
+      rc02  = rc0  * rc0   ;   rc02r = 1. / rc02
+      rc03  = rc02 * rc0
+      rc04  = rc03 * rc0
       rc03_sqrt2_galp = rc03 / SQRT(2._wp) / rn_clim_galp
-      rsbc_mb   = 0.5 * (15.8*rn_crban)**(2./3.) ! Surf. bound. cond. from Mellor and Blumberg
-      rsbc_std  = 3.75                           ! Surf. bound. cond. standard (prod=diss)
-      rsbc_tke1 = (-rsc_tke*rn_crban/(rcm_sf*ra_sf*rl_sf))**(2./3.) ! k_eps = 53.  Dirichlet + Wave breaking 
-      rsbc_tke2 = 0.5 / rau0
+      rsbc_mb   = 0.5_wp * (15.8_wp*rn_crban)**(2._wp/3._wp)               ! Surf. bound. cond. from Mellor and Blumberg
+      rsbc_std  = 3.75_wp                                                  ! Surf. bound. cond. standard (prod=diss)
+      rsbc_tke1 = (-rsc_tke*rn_crban/(rcm_sf*ra_sf*rl_sf))**(2._wp/3._wp)  ! k_eps = 53.  Dirichlet + Wave breaking 
+      rsbc_tke2 = 0.5_wp / rau0
       rsbc_tke3 = rdt * rn_crban                                                         ! Neumann + Wave breaking
-      rsbc_zs   = rn_charn/grav                                                          ! Charnock formula
+      rsbc_zs   = rn_charn / grav                                                        ! Charnock formula
       rsbc_psi1 = rc0**rpp * rsbc_tke1**rmm * rl_sf**rnn                           ! Dirichlet + Wave breaking
-      rsbc_psi2 = -0.5 * rdt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi                   ! Neumann + NO Wave breaking 
-      rsbc_psi3 = -0.5 * rdt * rc0**rpp * rl_sf**rnn / rsc_psi  * (rnn + rmm*ra_sf) ! Neumann + Wave breaking
-      rfact_tke = -0.5 / rsc_tke * rdt               ! Cst used for the Diffusion term of tke
-      rfact_psi = -0.5 / rsc_psi * rdt               ! Cst used for the Diffusion term of tke
-
-      ! Wall proximity function
+      rsbc_psi2 = -0.5_wp * rdt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi                   ! Neumann + NO Wave breaking 
+      rsbc_psi3 = -0.5_wp * rdt * rc0**rpp * rl_sf**rnn / rsc_psi  * (rnn + rmm*ra_sf) ! Neumann + 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
+      !                                !* set vertical eddy coef. to the background value
       DO jk = 1, jpk
          avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
@@ -1218,5 +1162,5 @@
          avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)
       END DO
-      !                               !* read or initialize all required files 
+      !                                !* read or initialize all required files 
       CALL gls_rst( nit000, 'READ' )
       !
@@ -1236,5 +1180,5 @@
      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
@@ -1327,10 +1271,8 @@
       WRITE(*,*) 'zdf_gls: You should not have seen this print! error?', kt
    END SUBROUTINE zdf_gls
-   SUBROUTINE zdf_gls_init          ! Empty routine
-   END SUBROUTINE zdf_gls_init
-   SUBROUTINE gls_rst( kt,cdrw )          ! Empty routine
+   SUBROUTINE gls_rst( kt, cdrw )          ! Empty routine
       INTEGER         , INTENT(in) ::   kt         ! ocean time-step
       CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
-      WRITE(*,*) 'gls_Rst: You should not have seen this print! error?', kt, cdrw
+      WRITE(*,*) 'gls_rst: You should not have seen this print! error?', kt, cdrw
    END SUBROUTINE gls_rst
 #endif
@@ -1338,2 +1280,3 @@
    !!======================================================================
 END MODULE zdfgls
+
