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 2293)
+++ branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90	(revision 2299)
@@ -59,52 +59,52 @@
    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
+
+   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) ::   ra1           =  0.92_wp
+   REAL(wp) ::   ra2           =  0.74_wp
+   REAL(wp) ::   rb1           = 16.60_wp
+   REAL(wp) ::   rb2           = 10.10_wp         
+   REAL(wp) ::   re2           =  1.33_wp         
+   REAL(wp) ::   rl1           =  0.107_wp
+   REAL(wp) ::   rl2           =  0.0032_wp
+   REAL(wp) ::   rl3           =  0.0864_wp
+   REAL(wp) ::   rl4           =  0.12_wp
+   REAL(wp) ::   rl5           = 11.9_wp
+   REAL(wp) ::   rl6           =  0.4_wp
+   REAL(wp) ::   rl7           =  0.0_wp
+   REAL(wp) ::   rl8           =  0.48_wp
+   REAL(wp) ::   rm1           =  0.127_wp
+   REAL(wp) ::   rm2           =  0.00336_wp
+   REAL(wp) ::   rm3           =  0.0906_wp
+   REAL(wp) ::   rm4           =  0.101_wp
+   REAL(wp) ::   rm5           = 11.2_wp
+   REAL(wp) ::   rm6           =  0.4_wp
+   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) ::   rnc_tke, rsc_psi, rpsi1, rpsi2, rpsi3, rsc_psi0
+   REAL(wp) ::   rpsi3m, rpsi3p, rpp, rmm, rnn
 
    !! * Substitutions
@@ -169,5 +169,5 @@
           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)
+          ustars2(ji,jj) = rsbc_tke2 * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)
           !
           ! bottom friction 
@@ -183,5 +183,5 @@
       ! Compute surface roughness length according to Charnock formula:
       IF (ln_crban) THEN
-        zhsro(:,:) = MAX(rn_sbc_zs * ustars2(:,:), hsro)
+        zhsro(:,:) = MAX(rsbc_zs * ustars2(:,:), hsro)
       ELSE
         zhsro(:,:) = hsro
@@ -200,5 +200,5 @@
                   &                            / (  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)
+               eps(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT(en(ji,jj,jk)) / mxln(ji,jj,jk)
             ENDDO
          ENDDO
@@ -218,5 +218,5 @@
                   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)
+                  zwall (ji,jj,jk) = ( 1. + re2 * ( zup / MAX( zdown, rsmall ) )**2. ) * tmask(ji,jj,jk)
                ENDDO
             ENDDO
@@ -256,11 +256,11 @@
                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
+               ! Compute a wall function from 1. to rsc_psi*zwall/rsc_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
+               ! 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) = rn_sc_psi / (zsigpsi * rn_sc_psi + (1.-zsigpsi) * rn_sc_psi0 / MAX(zwall(ji,jj,jk),1.))
+                 zwall_psi(ji,jj,jk) = rsc_psi / (zsigpsi * rsc_psi + (1.-zsigpsi) * rsc_psi0 / MAX(zwall(ji,jj,jk),1.))
                ELSE
                  zwall_psi(ji,jj,jk) = 1.
@@ -268,5 +268,5 @@
                !
                ! building the matrix
-               zcof = zfact_tke * tmask(ji,jj,jk)
+               zcof = rfact_tke * tmask(ji,jj,jk)
                !
                ! lower diagonal
@@ -294,5 +294,5 @@
         DO jj = 2, jpjm1
           DO ji = fs_2, fs_jpim1   ! vector opt.
-            zwall_psi(ji,jj,1) = rn_sc_psi/rn_sc_psi0
+            zwall_psi(ji,jj,1) = rsc_psi/rsc_psi0
           END DO
         END DO
@@ -310,5 +310,5 @@
       !                      ! balance between the production and the dissipation terms including the wave effect
       !
-      en(:,:,1) = MAX( rn_sbc_tke1 * ustars2(:,:), rn_emin )
+      en(:,:,1) = MAX( rsbc_tke1 * ustars2(:,:), rn_emin )
       z_elem_a(:,:,1) = en(:,:,1)
       z_elem_c(:,:,1) = 0.
@@ -316,5 +316,5 @@
       ! 
       ! one level below
-      en(:,:,2) = MAX( rn_sbc_tke1 * ustars2(:,:) * ( (zhsro(:,:)+fsdepw(:,:,2))/zhsro(:,:) )**rn_a_sf, rn_emin )
+      en(:,:,2) = MAX( rsbc_tke1 * ustars2(:,:) * ( (zhsro(:,:)+fsdepw(:,:,2))/zhsro(:,:) )**ra_sf, rn_emin )
       z_elem_a(:,:,2) = 0.
       z_elem_c(:,:,2) = 0.
@@ -325,5 +325,5 @@
       !                      ! balance between the production and the dissipation terms
       !
-      en(:,:,1) = MAX( rn_c02r * ustars2(:,:), rn_emin )
+      en(:,:,1) = MAX( rc02r * ustars2(:,:), rn_emin )
       z_elem_a(:,:,1) = en(:,:,1) 
       z_elem_c(:,:,1) = 0.
@@ -331,5 +331,5 @@
       !
       ! one level below
-      en(:,:,2) = MAX( rn_c02r * ustars2(:,:), rn_emin )
+      en(:,:,2) = MAX( rc02r * ustars2(:,:), rn_emin )
       z_elem_a(:,:,2) = 0.
       z_elem_c(:,:,2) = 0.
@@ -343,5 +343,5 @@
       !
       ! Dirichlet conditions at k=1 (Cosmetic)
-      en(:,:,1) = MAX( rn_sbc_tke1 * ustars2(:,:), rn_emin )
+      en(:,:,1) = MAX( rsbc_tke1 * ustars2(:,:), rn_emin )
       z_elem_a(:,:,1) = en(:,:,1)
       z_elem_c(:,:,1) = 0.
@@ -350,5 +350,5 @@
       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)
+      zflxs(:,:) = rsbc_tke3 * ustars2(:,:)**1.5 * ((zhsro(:,:)+fsdept(:,:,1))/zhsro(:,:) )**(1.5*ra_sf)
       en(:,:,2) = en(:,:,2) + zflxs(:,:)/fse3w(:,:,2)
       !
@@ -356,5 +356,5 @@
       !
       ! Dirichlet conditions at k=1 (Cosmetic)
-      en(:,:,1) = MAX( rn_c02r * ustars2(:,:), rn_emin )
+      en(:,:,1) = MAX( rc02r * ustars2(:,:), rn_emin )
       z_elem_a(:,:,1) = en(:,:,1)
       z_elem_c(:,:,1) = 0.
@@ -388,5 +388,5 @@
             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 )
+            en(ji,jj,ibot  ) = MAX( rc02r * ustarb2(ji,jj), rn_emin )
             !
             ! Just above last level, Dirichlet condition again
@@ -394,5 +394,5 @@
             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 ) 
+            en(ji,jj,ibotm1) = MAX( rc02r * ustarb2(ji,jj), rn_emin ) 
          END DO
       END DO
@@ -411,5 +411,5 @@
             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 )
+            en(ji,jj,ibot) = MAX( rc02r * ustarb2(ji,jj), rn_emin )
             !
             ! Just above last level: Neumann condition
@@ -482,5 +482,5 @@
          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))
+               psi(ji,jj,jk)  = SQRT(en(ji,jj,jk)) / (rc0 * mxln(ji,jj,jk))
             ENDDO
          ENDDO
@@ -492,5 +492,5 @@
          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 
+               psi(ji,jj,jk)  = rc02 * en(ji,jj,jk) * mxln(ji,jj,jk)**rnn 
             ENDDO
          ENDDO
@@ -516,14 +516,14 @@
                dir = 0.5 + sign(0.5,rn2(ji,jj,jk))
                !
-               rn_psi3 = dir*rn_psi3m+(1-dir)*rn_psi3p
+               rpsi3 = dir*rpsi3m+(1-dir)*rpsi3p
                !
                ! shear prod. - stratif. destruction
-               prod = rn_psi1 * zratio * shear(ji,jj,jk)
+               prod = rpsi1 * zratio * shear(ji,jj,jk)
                !
                ! stratif. destruction
-               buoy = rn_psi3 * 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 = rn_psi2 * zratio * zwall(ji,jj,jk) * eps(ji,jj,jk)
+               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)
@@ -533,5 +533,5 @@
                !                                                        
                ! building the matrix
-               zcof = zfact_psi * zwall_psi(ji,jj,jk) * tmask(ji,jj,jk)
+               zcof = rfact_psi * zwall_psi(ji,jj,jk) * tmask(ji,jj,jk)
                ! lower diagonal
                z_elem_a(ji,jj,jk) = zcof * ( avm  (ji,jj,jk  ) + avm  (ji,jj,jk-1) )   &
@@ -563,6 +563,6 @@
       !                      ! 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)
+      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.
@@ -570,7 +570,7 @@
       !
       ! 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)
+      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.
@@ -582,5 +582,5 @@
       !
       zdep(:,:) = vkarmn * zhsro(:,:)
-      psi (:,:,1) = rn_c0**zpp * en(:,:,1)**zmm * zdep(:,:)**znn * tmask(:,:,1)
+      psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
       z_elem_a(:,:,1) = psi(:,:,1)
       z_elem_c(:,:,1) = 0.
@@ -589,5 +589,5 @@
       ! one level below
       zdep(:,:) = vkarmn * ( zhsro(:,:) + fsdepw(:,:,2) )
-      psi (:,:,2) = rn_c0**zpp * en(:,:,1)**zmm * zdep(:,:)**znn * tmask(:,:,1)
+      psi (:,:,2) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
       z_elem_a(:,:,2) = 0.
       z_elem_c(:,:,2) = 0.
@@ -600,6 +600,6 @@
       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)
+      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.
@@ -611,7 +611,7 @@
       !
       ! 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         
+      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)
       !
@@ -619,5 +619,5 @@
       !
       zdep(:,:) = vkarmn * zhsro(:,:)
-      psi (:,:,1) = rn_c0**zpp * en(:,:,1)**zmm * zdep(:,:)**znn * tmask(:,:,1)
+      psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
       z_elem_a(:,:,1) = psi(:,:,1)
       z_elem_c(:,:,1) = 0.
@@ -630,5 +630,5 @@
       ! 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.)
+      zflxs(:,:) = rsbc_psi2 * ( avm(:,:,1) + avm(:,:,2) ) * en(:,:,1)**rmm * zdep**(rnn-1.)
       psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2)
       !     
@@ -653,5 +653,5 @@
             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
+            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.
@@ -660,5 +660,5 @@
             ! 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
+            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.
@@ -679,5 +679,5 @@
             ! 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
+            psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn
             !
             z_elem_a(ji,jj,ibot) = 0.
@@ -691,6 +691,6 @@
             ! 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.)
+            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
@@ -734,5 +734,5 @@
          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)
+              eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * en(ji,jj,jk) * SQRT(en(ji,jj,jk)) / psi(ji,jj,jk)
             ENDDO
          ENDDO
@@ -754,5 +754,5 @@
          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) 
+              eps(ji,jj,jk) = rc04 * en(ji,jj,jk) * psi(ji,jj,jk) 
             ENDDO
          ENDDO
@@ -764,5 +764,5 @@
          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)
+              eps(ji,jj,jk) = rc0**(3.+rpp/rnn) * en(ji,jj,jk)**(1.5+rmm/rnn) * psi(ji,jj,jk)**(-1./rnn)
             ENDDO
          ENDDO
@@ -778,5 +778,5 @@
                ! 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)
+               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 )
@@ -801,14 +801,13 @@
                ! 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 = MIN( gh, rgh0   )
+               gh = MAX( gh, rghmin )
                ! 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)
+               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) = rn_c_diff * sh * tmask(ji,jj,jk)
-               avmv(ji,jj,jk) = rn_c_diff * sm * tmask(ji,jj,jk)
+               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
@@ -824,22 +823,21 @@
                ! 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
+               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*rn_f6
-               gm = MIN ( (rn_d0 - rn_d1*gh + rn_d3*gh**2.)/(rn_d2-rn_d4*gh) , gm )
+               gm = gm*rf6
+               gm = MIN ( (rd0 - rd1*gh + rd3*gh**2.)/(rd2-rd4*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
+               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) = rn_c_diff * sh * tmask(ji,jj,jk)
-               avmv(ji,jj,jk) = rn_c_diff * sm * tmask(ji,jj,jk)
+               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
@@ -852,5 +850,5 @@
       DO jj = 2, jpjm1
          DO ji = fs_2, fs_jpim1   ! vector opt.
-            avmv(ji,jj,1) = rn_cm_sf / SQRT(2.)
+            avmv(ji,jj,1) = rcm_sf / SQRT(2.)
          END DO
       END DO
@@ -859,5 +857,5 @@
             ibot=mbathy(ji,jj)
             ibotm1=ibot-1
-            avmv(ji,jj,ibot) = rn_c0 / SQRT(2.)
+            avmv(ji,jj,ibot) = rc0 / SQRT(2.)
          END DO
       END DO
@@ -988,21 +986,21 @@
       !
       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
+      rpp       = 0.
+      rmm       = 1.
+      rnn       = 1.
+      rsc_tke = 1.96
+      rsc_psi = 1.96
+      rpsi1 = 0.9
+      rpsi3p = 1.
+      rpsi2  = 0.5
       !
          SELECT CASE ( nn_stab_func )
          !
          CASE( 0, 1 )               ! G88 or KC stability functions
-            rn_psi3m = 2.53
+            rpsi3m = 2.53
          CASE( 2 )                  ! Canuto A stability functions
-            rn_psi3m = 2.38
+            rpsi3m = 2.38
          CASE( 3 )                  ! Canuto B stability functions
-            rn_psi3m = 2.38         ! caution : constant not identified 
+            rpsi3m = 2.38         ! caution : constant not identified 
          END SELECT
       !
@@ -1010,21 +1008,21 @@
       !
       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
+      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
-           rn_psi3m = -0.52
+           rpsi3m = -0.52
         CASE( 2 )                  ! Canuto A stability functions
-           rn_psi3m = -0.629 
+           rpsi3m = -0.629 
         CASE( 3 )                  ! Canuto B stability functions
-           rn_psi3m = -0.566
+           rpsi3m = -0.566
         END SELECT
       !
@@ -1032,21 +1030,21 @@
       !
       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
+      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
-           rn_psi3m = -0.58
+           rpsi3m = -0.58
         CASE( 2 )                  ! Canuto A stability functions
-           rn_psi3m = -0.64
+           rpsi3m = -0.64
         CASE( 3 )                  ! Canuto B stability functions
-           rn_psi3m = -0.64        ! caution : constant not identified
+           rpsi3m = -0.64        ! caution : constant not identified
         END SELECT
       !
@@ -1054,21 +1052,21 @@
       !
       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
+      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
-           rn_psi3m = 0.1
+           rpsi3m = 0.1
         CASE( 2 )                  ! Canuto A stability functions
-           rn_psi3m = 0.05
+           rpsi3m = 0.05
         CASE( 3 )                  ! Canuto B stability functions
-           rn_psi3m = 0.05         ! caution : constant not identified
+           rpsi3m = 0.05         ! caution : constant not identified
         END SELECT
       !
@@ -1083,70 +1081,70 @@
       !
       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
+      rc2 = 0.
+      rc3 = 0.
+      rc_diff = 1.
+      rc0     = 0.5544
+      rcm_sf = 0.9884
+      rghmin = -0.28
+      rgh0   = 0.0233
+      rghcri = 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
+      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
       !
       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
+      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
       !
       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
+      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
       !
       END SELECT
@@ -1157,21 +1155,21 @@
       ! 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)                                 
+         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)                                 
       ELSE
-         rn_sc_psi0 = rn_sc_psi
+         rsc_psi0 = rsc_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.) &
+      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.) &
                &                                       )
 
@@ -1182,40 +1180,40 @@
          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,*) 'Parameter  m = ',rmm
+         WRITE(numout,*) 'Parameter  n = ',rnn
+         WRITE(numout,*) 'Parameter  p = ',rpp
+         WRITE(numout,*) 'rpsi1   = ',rpsi1
+         WRITE(numout,*) 'rpsi2   = ',rpsi2
+         WRITE(numout,*) 'rpsi3m  = ',rpsi3m
+         WRITE(numout,*) 'rpsi3p  = ',rpsi3p
+         WRITE(numout,*) 'rsc_tke = ',rsc_tke
+         WRITE(numout,*) 'rsc_psi = ',rsc_psi
+         WRITE(numout,*) 'rsc_psi0 = ',rsc_psi0
+         WRITE(numout,*) 'rc0     = ',rc0
          WRITE(numout,*)
          WRITE(numout,*) 'Shear free turbulence parameters:'
-         WRITE(numout,*) 'rn_cm_sf  = ',rn_cm_sf
-         WRITE(numout,*) 'rn_a_sf   = ',rn_a_sf
-         WRITE(numout,*) 'rn_l_sf   = ',rn_l_sf
+         WRITE(numout,*) 'rcm_sf  = ',rcm_sf
+         WRITE(numout,*) 'ra_sf   = ',ra_sf
+         WRITE(numout,*) 'rl_sf   = ',rl_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) / rn_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
+      rc02r = 1. / rc0**2.
+      rc02  = rc0**2._wp
+      rc03  = rc0**3._wp
+      rc04  = rc0**4._wp
+      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_tke3 = rdt * rn_crban                                                         ! Neumann + Wave breaking
+      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
