   !!======================================================================

  SUBROUTINE init_turbcoef(ect,psi)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE init_turbcoef  ***
   !&E
   !&E ** Purpose : initialize coefficients of turbulent parametrization
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by : main
   !&E
   !&E ** External calls : 
   !&E
   !&E ** Used ij-arrays : nzinit, kzinit
   !&E
   !&E ** Modified variables : nz, kz
   !&E
   !&E ** Reference : 
   !&E  - John C. Warner - Ocean Modelling 2005 : 
   !&E       Performance of 4 turbulence closure models implemented using a GLS method
   !&E  - Lars Umlauf, Hans Burchard - Continental Shelf Research 2005 : 
   !&E       Second-order turbulence closure models for geophysical boundary layers. A review of recent work
   !&E  - Lars Umlauf, Hans Burchard, Karsten Bolding : 2007
   !&E        GOTM Version 4.0     
   !&E       
   !&E
   !&E => modele (formulate with Generic Length Scale method)
   !&E  turb_2eq_option = 1 : modele k-kl MY2.5
   !&E                    2 : modele k-epsilon
   !&E                    3 : modele k-omega
   !&E                    4 : modele generique (Umlauf Burchard 2003)
   !&E  
   !&E => Fonctions de stabilite  
   !&E  OptionSF = 1 : Gibson Launder 1978
   !&E             2 : Mellor Yamada  1982  
   !&E             3 : Kantha Clayson 1994
   !&E             4 : Luyten & al.   1996
   !&E             5 : Canuto & al. A 2001
   !&E             6 : Canuto & al. B 2001
   !&E             7 : Cheng  & al.   2002
   !&E  
   !&E ** History :
   !&E       ! 2004-08-17
   !&E       ! 2007-09-10  (Pascal Douillet) initialisation nbouc=1
   !&E       ! 2009-12-12  (Thomas Duhaut) initialisation for 2 equations (k-kl,keps,komega)
   !&E       ! 2013-05-01  (S. Petton) Reformulation of stability functions using a unifying notation
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE parameters
   USE comvars2d,ONLY:h0
   USE comvars3d,ONLY:nz,kz,nzinit,kzinit,turb_2eq_option,turb_nbeq,OptionSF
   USE comturb  ,ONLY:sigmaect,sigmapsi,c1turb,c2turb,c3pturb,ectmin,    &
                      psimin,epsmin,fwall,cmu,c3mturb,pturb,mturb,nturb, &
                      eps,emin,alpha_n_min,                              &
                      sf_d0,sf_d1,sf_d2,sf_d3,sf_d4,sf_d5,               &
                      sf_n0,sf_n1,sf_n2,sf_nb0,sf_nb1,sf_nb2,            &
                      lim_am0,lim_am1,lim_am2,lim_am3,lim_am4,lim_am5,lim_am6

   IMPLICIT NONE

   !! * Arguments
   REAL(kind=rsh),DIMENSION(kmax,limin:limax,ljmin:ljmax),INTENT(out) :: ect,psi

   !! * Local declarations
   INTEGER :: i,j
   REAL(kind=rsh),DIMENSION(1:7) :: c1,c2,c3,c4,c5,c6,cb1,cb2,cb3,cb4,cb5,cbb,lim_an
   REAL(kind=rsh) :: a1,a2,a3,a4,a5,ab1,ab2,ab3,ab4,ab5,N,Nb

   !! * Executable part
   OptionSF=5
   
   ! Umlauf, Burchard : Table 1
   c1  = (/ 3.6000_rsh, 6.0000_rsh, 6.0000_rsh, 3.0000_rsh, 5.0000_rsh, 5.0000_rsh, 5.0000_rsh /) 
   c2  = (/ 0.8000_rsh, 0.3200_rsh, 0.3200_rsh, 0.8000_rsh, 0.8000_rsh, 0.6983_rsh, 0.7983_rsh /)
   c3  = (/ 1.2000_rsh, 0.0000_rsh, 0.0000_rsh, 2.0000_rsh, 1.9680_rsh, 1.9664_rsh, 1.9680_rsh /)
   c4  = (/ 1.2000_rsh, 0.0000_rsh, 0.0000_rsh, 1.1180_rsh, 1.1360_rsh, 1.0940_rsh, 1.1360_rsh /)
   c5  = (/ 0.0000_rsh, 0.0000_rsh, 0.0000_rsh, 0.0000_rsh, 0.0000_rsh, 0.0000_rsh, 0.0000_rsh /)
   c6  = (/ 0.5000_rsh, 0.0000_rsh, 0.0000_rsh, 0.5000_rsh, 0.4000_rsh, 0.4950_rsh, 0.5000_rsh /)
   cb1 = (/ 3.0000_rsh, 3.7280_rsh, 3.7280_rsh, 3.0000_rsh, 5.9500_rsh, 5.6000_rsh, 5.5200_rsh /)
   cb2 = (/ 0.3333_rsh, 0.0000_rsh, 0.7000_rsh, 0.3333_rsh, 0.6000_rsh, 0.6000_rsh, 0.2134_rsh /)
   cb3 = (/ 0.3333_rsh, 0.0000_rsh, 0.7000_rsh, 0.3333_rsh, 1.0000_rsh, 1.0000_rsh, 0.3570_rsh /)
   cb4 = (/ 0.0000_rsh, 0.0000_rsh, 0.0000_rsh, 0.0000_rsh, 0.0000_rsh, 0.0000_rsh, 0.0000_rsh /)
   cb5 = (/ 0.3333_rsh, 0.0000_rsh, 0.2000_rsh, 0.3333_rsh, 0.3333_rsh, 0.3333_rsh, 0.3333_rsh /)
   cbb = (/ 0.8000_rsh, 0.6102_rsh, 0.6102_rsh, 0.8000_rsh, 0.7200_rsh, 0.4770_rsh, 0.8200_rsh /)

   ! Umlauf, Burchard : Equations (A.10) et (A.12)
   a1  = 2.0_rsh/3.0_rsh - c2(OptionSF)/2.0_rsh
   a2  = 1.0_rsh - c3(OptionSF)/2.0_rsh
   a3  = 1.0_rsh - c4(OptionSF)/2.0_rsh
   a4  = c5(OptionSF)/2.0_rsh
   a5  = 1.0_rsh/2.0_rsh - c6(OptionSF)/2.0_rsh
   ab1 = 1.0_rsh - cb2(OptionSF)
   ab2 = 1.0_rsh - cb3(OptionSF)
   ab3 = 2.0_rsh * ( 1.0_rsh - cb4(OptionSF))
   ab4 = 2.0_rsh * ( 1.0_rsh - cb5(OptionSF))
   ab5 = 2.0_rsh*cbb(OptionSF)*( 1.0_rsh - cb5(OptionSF))
   N   = c1(OptionSF)/2.0_rsh
   Nb  = cb1(OptionSF)

   ! Umlauf, Burchard : Equations (A.9)
   sf_d0  =   36.0_rsh* N**3.0_rsh * Nb**2.0_rsh
   sf_d1  =   84.0_rsh*a5*ab3* N**2.0_rsh*Nb + 36.0_rsh*ab5* N**3.0_rsh*Nb
   sf_d2  =   9.0_rsh*(ab2**2.0_rsh-ab1**2.0_rsh)* N**3.0_rsh - 12.0_rsh*(a2**2.0_rsh-3.0_rsh*a3**2.0_rsh)*N*Nb**2.0_rsh
   sf_d3  =   12.0_rsh*a5*ab3*(a2*ab1-3.0_rsh*a3*ab2)* N      &
            + 12.0_rsh*a5*ab3*(a3**2.0_rsh-a2**2.0_rsh)* Nb   &
            + 12.0_rsh*ab5*(3.0_rsh*a3**2.0_rsh-a2**2.0_rsh)*N*Nb
   sf_d4  =   48.0_rsh*a5**2.0_rsh*ab3**2.0_rsh*N + 36.0_rsh*a5*ab3*ab5* N**2.0_rsh
   sf_d5  =   3.0_rsh*(a2**2.0_rsh-3.0_rsh*a3**2.0_rsh)*(ab1**2.0_rsh-ab2**2.0_rsh)*N


   sf_n0  =   36.0_rsh*a1* N**2.0_rsh* Nb**2.0_rsh
   sf_n1  = - 12.0_rsh*a5*ab3*(ab1+ab2)* N**2.0_rsh           &
            + 8.0_rsh*a5*ab3*(6.0_rsh*a1-a2-3.0_rsh*a3)*N*Nb  &
            + 36.0_rsh*a1*ab5* N**2.0_rsh*Nb
   sf_n2  =   9.0_rsh*a1*(ab2**2.0_rsh-ab1**2.0_rsh) * N**2.0_rsh
   sf_nb0 =   12.0_rsh*ab3 * N**3.0_rsh * Nb
   sf_nb1 =   12.0_rsh*a5*ab3**2.0_rsh  * N**2.0_rsh
   sf_nb2 =   9.0_rsh*a1*ab3*(ab1-ab2) * N**2.0_rsh + (  6.0_rsh*a1*(a2-3.0_rsh*a3) &
            - 4.0_rsh*(a2**2.0_rsh-3.0_rsh*a3**2.0_rsh) )*ab3 * N * Nb
   
   ! Umlauf, Burchard : Equation (47) et Table 3
   lim_an=(/-1.044_rsh, -1.120_rsh, -1.312_rsh, -0.993_rsh, -2.324_rsh, -3.093_rsh, -2.029_rsh /)
   alpha_n_min = lim_an(OptionSF)

   ! Umlauf, Burchard : Equations (43) et (44)
   lim_am0 = sf_d0*sf_n0 
   lim_am1 = sf_d0*sf_n1 + sf_d1*sf_n0
   lim_am2 = sf_d1*sf_n1 + sf_d4*sf_n0
   lim_am3 = sf_d4*sf_n1
   lim_am4 = sf_d2*sf_n0
   lim_am5 = sf_d2*sf_n1+sf_d3*sf_n0
   lim_am6 = sf_d3*sf_n1

   ! GOTM Version 4.0 : Equation (122)
   a1    =  2.0_rsh/3.0_rsh - c2(OptionSF)/2.0_rsh
   a3    =  1.    - c4(OptionSF)/2.0_rsh
   N     =  c1(OptionSF)/2.0_rsh
   cmu   =  ( (a2**2.0_rsh - 3.0_rsh*a3**2.0_rsh + 3.0_rsh*a1*N)/(3.0_rsh* N**2.0_rsh) )**0.25_rsh
   ! RM : Dans GOTM : c3mturb = compute_cpsi3.F90
   
   ! John C. Warner : Tables 1 & 2
   IF(turb_2eq_option==1) THEN ! modele k-kl MY2.5
     pturb=0.0_rsh
     mturb=1.0_rsh
     nturb=1.0_rsh
     sigmaect=2.44_rsh
     sigmapsi=2.44_rsh
     c1turb=0.90_rsh
     c2turb=0.50_rsh
     c3pturb=1.00_rsh
     IF (OptionSF == 3) THEN
       c3mturb=2.53_rsh        
     ELSEIF (OptionSF == 5) THEN
       c3mturb=2.38_rsh
     ELSE 
       c3mturb=2.53_rsh        
     END IF
     ectmin=1.0e-8
     psimin=1.0e-8
     epsmin=1.0e-12
     fwall=1.0_rsh
   END IF

   IF(turb_2eq_option==2) THEN ! modele k-epsilon
     pturb= 3.0_rsh
     mturb= 1.5_rsh
     nturb=-1.0_rsh
     sigmaect=1.0_rsh
     sigmapsi=1.3_rsh
     c1turb=1.44_rsh
     c2turb=1.92_rsh
     c3pturb=1.0_rsh
     IF (OptionSF == 3) THEN
       c3mturb=-0.52_rsh        
     ELSEIF (OptionSF == 5) THEN
       c3mturb=-0.63_rsh
     ELSEIF (OptionSF == 6) THEN
       c3mturb=-0.57_rsh
     ELSE 
       c3mturb=-0.52_rsh        
     END IF
     ectmin=1.0e-8 
     psimin=1.0e-12   
     epsmin=1.0e-12
     fwall=1.0_rsh
   END IF

   IF(turb_2eq_option==3) THEN ! modele k-omega
     pturb=-1.0_rsh
     mturb= 0.5_rsh
     nturb=-1.0_rsh
     sigmaect=2.0_rsh
     sigmapsi=2.0_rsh
     c1turb=0.555_rsh 
     c2turb=0.833_rsh 
     c3pturb=1.0_rsh
     IF (OptionSF == 3) THEN
       c3mturb=-0.58_rsh        
     ELSEIF (OptionSF == 5) THEN
       c3mturb=-0.64_rsh
     ELSE
       c3mturb=-0.58_rsh        
     END IF
     ectmin=1.0e-8 
     psimin=1.0e-12
     epsmin=1.0e-12 ! a modifier
     fwall=1.0_rsh
   END IF

  IF(turb_2eq_option==4) THEN ! modele Gen 
     pturb= 2.0_rsh
     mturb= 1.0_rsh
     nturb=-0.66666667_rsh
     sigmaect=0.8_rsh 
     sigmapsi=1.07_rsh
     c1turb=1.00_rsh
     c2turb=1.22_rsh
     c3pturb=1.00_rsh
     IF (OptionSF == 3) THEN
       c3mturb=0.10_rsh        
     ELSEIF (OptionSF == 5) THEN
       c3mturb=0.05_rsh
     ELSE
       c3mturb=0.10_rsh        
     END IF
     ectmin=1.0e-8 
     psimin=1.0e-12
     epsmin=1.0e-12 ! a modifier
     fwall=1.0_rsh
   END IF

!$OMP parallel
!$OMP do private(j) 
   DO j=ljmin,ljmaxp1
     DO i=limin,limaxp1
       IF (h0(i,j)/=-valmanq) THEN
         nz(:,i,j)=nzinit
       ELSE
         nz(:,i,j)=valmanq
       END IF
     ENDDO
   ENDDO
!$OMP end do

!$OMP do private(j) 
   DO j=ljminm1,ljmaxp1
     DO i=liminm1,limaxp1
       IF (h0(i,j)/=-valmanq) THEN
         kz(:,i,j)=kzinit
       ELSE
         kz(:,i,j)=valmanq
       END IF
     ENDDO
   ENDDO
!$OMP end do

!$OMP do private(j) 
   DO j=ljmin,ljmax
     DO i=limin,limax
       IF (h0(i,j)/=-valmanq) THEN
         IF (turb_nbeq==1) THEN
           emin=1.0e-8
           ect(:,i,j)=emin
         ELSE
           ect(:,i,j)=ectmin
         ENDIF
         psi(:,i,j)=psimin
         eps(:,i,j)=epsmin
       ELSE
         ect(:,i,j)=valmanq
         psi(:,i,j)=valmanq
         eps(:,i,j)=valmanq
       END IF
     ENDDO
   ENDDO
!$OMP end do
!$OMP end parallel

 END SUBROUTINE init_turbcoef
