
  SUBROUTINE turbpsi(xeuv,uz,vz,pn2,psi,ect)


!&E----------------------------------------------------------------------------------------------------
!&E                 ***  ROUTINE turbpsi  ***
!&E
!&E ** Purpose : compute the vertical viscosity and diffusion coefficients,
!&E              using the k parameterisation
!&E
!&E ** Description : compute the turbulent kinetic energy (ect),
!&E                  estimate nz and kz from ect and from an anlytical
!&E                  formulation of the mixing length (l)
!&E
!&E ** Called by : stepu, stepv
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays : h0,xeuv,z0b,nz,uz,vz,eps,ect,pn2,psi,tosx,tosy
!&E
!&E ** Modified variables : psi,nz
!&E
!&E ** History :
!&E       !  2004-09-15
!&E       !  2008-12    (V. Garnier) shift of sigma values : -1 < sig < 0
!&E       !  2009-03    (T. Duhaut) corrections
!&E       !  2009-04    (A. Mano, V. Garnier) MPI2D
!&E       !  2009-06    (T. Duhaut) Choice of BC (surface and bottom): Flux or Dirichlet (default)
!&E       !  2009-06    (F. Dumas) generalized sigma coordinates
!&E       !  2009-12    (T. Duhaut, F. Dumas) evolution
!&E       !  2011-01    (A.-C. Bennis) correction dtm instead of dt
!&E       !  2011-10    (V. Garnier) OMP
!&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
!&E       !  2011-12  (S. Petton, F. Dumas) Bottom stress : removing of min, max
!&E       !  2012-02  (A.C Bennis, F. Dumas) Mixing length due to surface waves and bottom rugosity
!&E                   Warner et al.  (2005) Ocean Modelling : Performance of four turbulence closure models...
!&E       !  2012-02  (F. Dumas) Initialization of nz, kz at kmax (+ other corrections)
!&E       !  2013-08  (A-C. Bennis and F. Dumas) wave-current coupling
!&E       !  2014-07  (F. Ganthy) Introduction of obstruction module
!&E       !  2017-09  (F. Ganthy) Turbulence advection (k-epsilon)
!&E       !  2017-10  (F. Ganthy) Minor change for obstructions (remove call to
!&E       !                       obstruction-induced turbulence coefficient)
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   USE parameters
   USE comvars2d
   USE comvars3d


   USE cometdyn
   USE comturb


   USE toolmpi,  ONLY : ex_i_rsh,ex_j_rsh





   IMPLICIT NONE

!! * Arguments
   REAL(kind=rsh),DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2),INTENT(in)    :: xeuv

   REAL(kind=rsh),DIMENSION(kmax,liminm2:limaxp2,ljminm2:ljmaxp2),INTENT(in)    :: uz,vz
   REAL(kind=rsh),DIMENSION(kmax,limin  :limax  ,ljmin  :ljmax  ),INTENT(in)    :: pn2
   REAL(kind=rsh),DIMENSION(kmax,limin  :limax  ,ljmin  :ljmax  ),INTENT(inout) :: psi,ect

!! * Local declarations

   INTEGER        :: k,i,j,kkmax

   REAL(kind=rsh) :: nzms,nzmi,kfond
   REAL(kind=rsh) :: sqrtcmu,dtm,hp
   REAL(kind=rsh) :: cd,prod,dest,enplus
   REAL(kind=rsh),DIMENSION(kmax)   :: a,bb,c,y,u0,p0
   REAL(kind=rsh),DIMENSION(0:kmax) :: q0

   REAL(kind=rsh) :: kZb,kZs,FpsiS,db,ds


!!----------------------------------------------------------------------
!! * Executable part

!$OMP PARALLEL DEFAULT(none) SHARED(psi,ect,eps,h0,xeuv,uz,vz,nz,kz,pn2,    &
!$OMP& limin,limax,ljminp2,ljmaxm1,ig,id,dsigu,dsigw,                       &
!$OMP& z0b,z0s,turb_2eq_option,cmu,hm,halfdt,sigmapsi,psimin,icon,          &





!$OMP& sig,sigw,                                                            &

!$OMP& c1turb,c2turb,c3mturb,c3pturb)                                       &
!$OMP& PRIVATE(i,j,k,gamma,l_bc_turb,kkmax,sqrtcmu,dtm,hp,enplus,           &
!$OMP& kZb,kZs,cd,kfond,db,ds,fwall,nzms,nzmi,dest,prod,FpsiS,              &



!$OMP& a,bb,c,y,q0,p0,u0)

!------------------------------------------------------------------------------
!  coefficients de dispersion horizontale (pour calcul d advection dispersion)
!------------------------------------------------------------------------------
!
! dans l ideal il faut gerer c1, c2 et c3 et gamma devra disparaitre
   IF ( turb_2eq_option == 1) gamma = 1.0_rsh ! k-kl
   IF ( turb_2eq_option == 2) gamma = 1.0_rsh ! k-epsilon
   IF ( turb_2eq_option == 3) gamma = 1.0_rsh ! k-omega
   IF ( turb_2eq_option == 4) gamma = 1.0_rsh ! generic length scale

   l_bc_turb=.true.       ! Von Neumann boundary condition (we prescribe a flux in the middle of the surface or bottom cell)
! l_bc_turb=.false.    ! Dirichlet boundary condition (we prescribe the value at the interface (bottom or surface)

   IF (l_bc_turb) THEN
     kkmax=kmax-1
   ELSE
     kkmax=kmax
   END IF

   sqrtcmu =SQRT(cmu)
   fwall=1.0_rsh      ! private variable, so must be initialized

   dtm=halfdt

!$OMP DO SCHEDULE(RUNTIME)
   DO j=ljminp2,ljmaxm1
     DO i=MAX0(limin,ig(j)+1),MIN0(limax,id(j)-1)

       hp=h0(i,j)+xeuv(i,j)
       IF (hp > hm) THEN


! What a mess !!!! grrrr ! cval
! for the surface : Dirichlet BC uses kZb nd Neumann BC uses FpsiS
! for the bottom : Dirichlet and Neumann BC use the same variable kZb

! kappa*z at the bottom
         IF ( l_bc_turb) THEN          ! Neumann bc (flux) at u position

           kZb = 0.4_rsh*((sig(1)+1.0_rsh)*hp+z0b(i,j))

         ELSE                          ! Dirichlet (value) at the interface (bottom)
           kZb = 0.4_rsh*z0b(i,j)
         END IF

! kappa*z at the surface
         kZs = 0.4_rsh*z0s(i,j)

! nu/sigmapsi
         nz(:,i,j)=nz(:,i,j)/sigmapsi

         k=1
!*****************************


! (sig(1)+1.0_rsh)*hp est la distance depuis le fond
!    cd=min((0.4_rsh/log(((sig(1)+1.0_rsh)*hp+z0b(i,j))/z0b(i,j)))**2,0.1_rsh)
!cd=MIN((0.4_rsh/LOG((sig(1)+1.0_rsh)*hp/z0b(i,j)))**2,0.1_rsh)
         cd=(0.4_rsh/LOG((sig(k)+1.0_rsh)*hp/z0b(i,j)+1.0_rsh))**2





! kfond = epsilon at the bottom (for Dirichlet BC at least and for Neumann BC d epsilon /dz)....
! No comment for the choice of the name.... ! cval
         IF ( turb_2eq_option == 1 ) THEN   ! k-kl
           IF ( l_bc_turb) THEN
             kfond = cd*(0.5_rsh*(uz(k,i,j)**2+uz(k,i-1,j)**2)               &
                        +0.5_rsh*(vz(k,i,j)**2+vz(k,i,j-1)**2))*0.4_rsh/(cmu**2)
           ELSE
             kfond = cd*(0.5_rsh*(uz(k,i,j)**2+uz(k,i-1,j)**2)               &
                        +0.5_rsh*(vz(k,i,j)**2+vz(k,i,j-1)**2))*kZb/(cmu**2)
           END IF

           db=hp*dsigu(k)

           ds=hp-db
           fwall = 1.0_rsh+1.33_rsh*( psi(k,i,j)*(hp/(db*ds))  &
                                /(ect(k,i,j)*0.4_rsh) )**2

         ELSE IF ( turb_2eq_option == 2 ) THEN   ! k-epsilon
           IF ( l_bc_turb) THEN
             kfond =-0.4_rsh*( (cd*(0.5_rsh*(uz(k,i,j)**2+uz(k,i-1,j)**2)                      &
                                   +0.5_rsh*(vz(k,i,j)**2+vz(k,i,j-1)**2)))**(1.5_rsh) )/kZb**2
           ELSE
             kfond = ( (cd*(0.5_rsh*(uz(k,i,j)**2+uz(k,i-1,j)**2)                              &
                           +0.5_rsh*(vz(k,i,j)**2+vz(k,i,j-1)**2)))**(1.5_rsh) )/kZb
           END IF
         ELSE IF ( turb_2eq_option == 3 ) THEN   ! k-omega
           IF ( l_bc_turb) THEN
             kfond =-0.4_rsh*SQRT( cd*(0.5_rsh*(uz(k,i,j)**2+uz(k,i-1,j)**2)                   &
                                      +0.5_rsh*(vz(k,i,j)**2+vz(k,i,j-1)**2)) )/kZb**2/(cmu**2)
           ELSE
             kfond = SQRT( cd*(0.5_rsh*(uz(k,i,j)**2+uz(k,i-1,j)**2)                           &
                              +0.5_rsh*(vz(k,i,j)**2+vz(k,i,j-1)**2)) )/kZb/(cmu**2)
           END IF
         ELSE IF ( turb_2eq_option == 4 ) THEN    ! generic length scale
           IF ( l_bc_turb) THEN
             kfond =-0.4_rsh*(2.0_rsh/3.0_rsh)*( cd*(0.5_rsh*(uz(k,i,j)**2+uz(k,i-1,j)**2)    &
                                                    +0.5_rsh*(vz(k,i,j)**2+vz(k,i,j-1)**2)) ) &
                                              /kZb**(5.0_rsh/3.0_rsh)
           ELSE
             kfond = SQRT( cd*(0.5_rsh*(uz(k,i,j)**2+uz(k,i-1,j)**2)                          &
                              +0.5_rsh*(vz(k,i,j)**2+vz(k,i,j-1)**2)) )/kZb**(2.0_rsh/3.0_rsh)
           END IF
         END IF

         nzms=(dsigw(k+1)*nz(k,i,j)+dsigw(k)*nz(k+1,i,j))/(2.0_rsh*dsigu(k+1))
         nzmi=nz(k,i,j)




         c(k)=-dtm*nzms/(hp*hp*dsigu(k+1))

         a(k)=0.0_rsh




         prod=c1turb*sigmapsi*nz(k,i,j)*0.5_rsh*((vz(k+1,i,j  )-vz(k,i,j  ))**2   &
                                +(vz(k+1,i,j-1)-vz(k,i,j-1))**2                   &
                                +(uz(k+1,i,j  )-uz(k,i,j  ))**2                   &
                                +(uz(k+1,i-1,j)-uz(k,i-1,j))**2)/(hp*dsigw(k))**2

         dest=0.0_rsh
         IF (icon /=0) THEN
           IF (pn2(k,i,j) < 0.0_rsh) dest=c3pturb*kz(k,i,j)*pn2(k,i,j)
           IF (pn2(k,i,j) > 0.0_rsh) dest=c3mturb*kz(k,i,j)*pn2(k,i,j)
         END IF



         IF (l_bc_turb) THEN

           enplus=dtm*nzmi*kfond/hp

           bb(k)=dsigw(k)-c(k)+dtm*dsigw(k)*c2turb*eps(k,i,j)*gamma*fwall/ect(k,i,j)
           IF (psi(k,i,j)<1.e-08) psi(k,i,j)=1e-08
           IF (ect(k,i,j)<1.e-08) ect(k,i,j)=1e-08
           y(k)=psi(k,i,j)*dsigw(k)+dtm*dsigw(k)*(prod-dest)*psi(k,i,j)/ect(k,i,j) &
                -enplus+dtm*dsigw(k)*c2turb*eps(k,i,j)*fwall*psi(k,i,j)/ect(k,i,j)*(gamma-1.0_rsh) 
         ELSE

           bb(k)=dsigw(k)-c(k)+dtm*nzmi/(hp*hp*dsigu(k)*0.5_rsh)        &
                 +dtm*dsigw(k)*c2turb*eps(k,i,j)*gamma*fwall/ect(k,i,j)
           enplus=dtm*nzmi*kfond/(hp*hp*dsigu(k)*0.5_rsh)

           y(k)=psi(k,i,j)*dsigw(k)+dtm*dsigw(k)*(prod-dest)*psi(k,i,j)/ect(k,i,j)                &
                +enplus+dtm*dsigw(k)*c2turb*eps(k,i,j)*fwall*psi(k,i,j)/ect(k,i,j)*(gamma-1.0_rsh) 
         END IF






         DO k=2,kkmax-1
!*****************************

           IF ( turb_2eq_option == 1 ) THEN

             db=db+hp*dsigu(k)

             ds=hp-db
             fwall = 1.0_rsh+1.33_rsh*( psi(k,i,j)*(hp/(db*ds))  &
                                   /(ect(k,i,j)*0.4_rsh) )**2
           END IF

           nzms=(dsigw(k+1)*nz(k,i,j)+dsigw(k)*nz(k+1,i,j))/(2.0_rsh*dsigu(k+1))




           c(k)=-dtm*nzms/(hp*hp*dsigu(k+1))




           a(k)=c(k-1)
           bb(k)=dsigw(k)-a(k)-c(k)+dtm*dsigw(k)*c2turb*eps(k,i,j)*gamma*fwall/ect(k,i,j)

           prod=c1turb*sigmapsi*nz(k,i,j)*0.5_rsh*((vz(k+1,i,j  )-vz(k,i,j  ))**2                   &
                                                  +(vz(k+1,i,j-1)-vz(k,i,j-1))**2                   &
                                                  +(uz(k+1,i,j  )-uz(k,i,j  ))**2                   &
                                                  +(uz(k+1,i-1,j)-uz(k,i-1,j))**2)/(hp*dsigw(k))**2

           dest=0.0_rsh
           IF (icon /=0) THEN
             IF (pn2(k,i,j) < 0.0_rsh) dest=c3pturb*kz(k,i,j)*pn2(k,i,j)
             IF (pn2(k,i,j) > 0.0_rsh) dest=c3mturb*kz(k,i,j)*pn2(k,i,j)
           END IF

           y(k)=psi(k,i,j)*dsigw(k)+dtm*dsigw(k)*(prod-dest)*psi(k,i,j)/ect(k,i,j)     &
                       +dtm*dsigw(k)*c2turb*eps(k,i,j)*fwall*psi(k,i,j)/ect(k,i,j)*(gamma-1.0_rsh)  






         END DO

!******************
         k=kkmax
!******************

         IF ( l_bc_turb) THEN

           prod=c1turb*sigmapsi*nz(k,i,j)*0.5_rsh*((vz(k+1,i,j  )-vz(k,i,j  ))**2                   &
                                                  +(vz(k+1,i,j-1)-vz(k,i,j-1))**2                   &
                                                  +(uz(k+1,i,j  )-uz(k,i,j  ))**2                   &
                                                  +(uz(k+1,i-1,j)-uz(k,i-1,j))**2)/(hp*dsigw(k))**2

           dest=0.0_rsh
           IF (icon /=0) THEN
             IF (pn2(k,i,j) < 0.0_rsh) dest=c3pturb*kz(k,i,j)*pn2(k,i,j)
             IF (pn2(k,i,j) > 0.0_rsh) dest=c3mturb*kz(k,i,j)*pn2(k,i,j)
           END IF


           nzms=(dsigw(k+1)*nz(k,i,j)+dsigw(k)*nz(k+1,i,j))/(2.0_rsh*dsigu(k+1))



           IF ( turb_2eq_option == 1 ) THEN

             db=db+hp*dsigu(k)

             ds=hp-db
             fwall = 1.0_rsh+1.33_rsh*( psi(k,i,j)*(hp/(db*ds))  &
                                   /(ect(k,i,j)*0.4_rsh) )**2

             FpsiS = -ect(kmax,i,j)*0.4_rsh
           ELSE IF ( turb_2eq_option == 2 ) THEN

!             FpsiS = (cmu**3)*(ect(kmax,i,j)**1.5_rsh)/0.4_rsh/((hp*sigw(kmax-1)*0.5_rsh)**2)
             FpsiS = (cmu**3)*(ect(kmax,i,j)**1.5_rsh)/0.4_rsh/((-1.0_rsh*hp*sigw(kmax-1)*0.5_rsh+z0s(i,j))**2)


           ELSE IF ( turb_2eq_option == 3 ) THEN
!cval pourquoi n y a t il pas waveqeuler la ?

             FpsiS = SQRT(ect(kmax,i,j))/0.4_rsh/((-hp*sigw(kmax-1)*0.5_rsh+z0s(i,j))**2)/cmu

           ELSE IF ( turb_2eq_option == 4 ) THEN
!cval pourquoi n y a t il pas waveqeuler la ?

             FpsiS = (2.0_rsh/3.0_rsh)*(cmu**2)*ect(kmax,i,j)/(0.4_rsh**(2.0_rsh/3.0_rsh)) &
                    /((-hp*sigw(kmax-1)*0.5_rsh+z0s(i,j))**(5.0_rsh/3.0_rsh))

           END IF

           c(k)=0.0_rsh
           a(k)=c(k-1)
           bb(k)=dsigw(k)-a(k)+dtm*dsigw(k)*c2turb*eps(k,i,j)*fwall/ect(k,i,j)*gamma


           y(k)=psi(k,i,j)*dsigw(k)+dtm*dsigw(k)*(prod-dest)*psi(k,i,j)/ect(k,i,j)         &
                +dtm*nzms*FpsiS/hp                                                         &
                +dtm*dsigw(k)*c2turb*eps(k,i,j)*fwall*psi(k,i,j)/ect(k,i,j)*(gamma-1.0_rsh) 







!*******************
           k=kmax
!*******************

           a(k)=0.0_rsh
           bb(k)=1.0_rsh
           IF (turb_2eq_option==1) THEN
             y(k) = ect(k,i,j)*kZs
           ELSE IF (turb_2eq_option==2) THEN

             y(k) = (cmu**3)*(ect(k,i,j)**1.5_rsh)/kZs




           ELSE IF (turb_2eq_option==3) THEN
             y(k) = SQRT(ect(k,i,j))/(kZs*cmu)
           ELSE IF (turb_2eq_option==4) THEN
             y(k) = (cmu**2)*ect(k,i,j)/(kZs**(2.0_rsh/3.0_rsh))
           END IF

         ELSE
! same code than above since kkmax=kmax if we use dirichlet boundary condition (l_bc_turb=.false.)
           a(k)=0.0_rsh
           bb(k)=1.0_rsh
           IF (turb_2eq_option==1) THEN
             y(k) = ect(k,i,j)*kZs
           ELSE IF (turb_2eq_option==2) THEN
             y(k) = (cmu**3)*(ect(k,i,j)**1.5_rsh)/kZs
           ELSE IF (turb_2eq_option==3) THEN
             y(k) = SQRT(ect(k,i,j))/(kZs*cmu)
           ELSE IF (turb_2eq_option==4) THEN
             y(k) = (cmu**2)*ect(k,i,j)/(kZs**(2.0_rsh/3.0_rsh))
           END IF
         END IF

!********************************
!  inversion du systeme 1DV
!********************************

!******************************
!  factorisation LU
!******************************
         q0(1)=-c(1)/bb(1)
         DO k=2,kmax-1
           p0(k)=a(k)*q0(k-1)+bb(k)
           q0(k)=-c(k)/p0(k)
         END DO
         p0(kmax)=a(kmax)*q0(kmax-1)+bb(kmax)

!******************************
!  resolution Lu0=y
!******************************
         u0(1)=y(1)/bb(1)
         DO k=2,kmax
           u0(k)=(y(k)-a(k)*u0(k-1))/p0(k)
         END DO

!******************************
!  resolution Upsi=u0
!******************************
         psi(kmax,i,j)=u0(kmax)
         DO k=kmax-1,1,-1
           psi(k,i,j)=q0(k)*psi(k+1,i,j)+u0(k)
         END DO



         psi(:,i,j)=MAX(psi(:,i,j),psimin)
         nz(:,i,j)=nz(:,i,j)*sigmapsi       ! since we have divided by sigmapsi at the entrance of the routine
       END IF
     END DO
   END DO
!$OMP END DO

!$OMP END PARALLEL


!-----------------------------------------------------------------------------

  END SUBROUTINE turbpsi
