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

! This file is part of Code_Saturne, a general-purpose CFD tool.
!
! Copyright (C) 1998-2021 EDF S.A.
!
! This program is free software; you can redistribute it and/or modify it under
! the terms of the GNU General Public License as published by the Free Software
! Foundation; either version 2 of the License, or (at your option) any later
! version.
!
! This program is distributed in the hope that it will be useful, but WITHOUT
! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
! FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
! details.
!
! You should have received a copy of the GNU General Public License along with
! this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
! Street, Fifth Floor, Boston, MA 02110-1301, USA.

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

!===============================================================================
! Function :
! --------

!> \file clsyvt.f90
!>
!> \brief Symmetry boundary conditions for vectors and tensors.
!>
!> Correspond to the code icodcl(ivar) = 4.
!>
!> Please refer to the
!> <a href="../../theory.pdf#clsyvt"><b>clsyvt</b></a> section of the
!> theory guide for more informations.
!-------------------------------------------------------------------------------

!-------------------------------------------------------------------------------
! Arguments
!______________________________________________________________________________.
!  mode           name          role                                           !
!______________________________________________________________________________!
!> \param[in]     nscal         total number of scalars
!> \param[in,out] icodcl        face boundary condition code:
!>                               - 1 Dirichlet
!>                               - 3 Neumann
!>                               - 4 sliding and
!>                                 \f$ \vect{u} \cdot \vect{n} = 0 \f$
!>                               - 5 smooth wall and
!>                                 \f$ \vect{u} \cdot \vect{n} = 0 \f$
!>                               - 6 rough wall and
!>                                 \f$ \vect{u} \cdot \vect{n} = 0 \f$
!>                               - 9 free inlet/outlet
!>                                 (input mass flux blocked to 0)
!> \param[in,out] rcodcl        boundary condition values:
!>                               - rcodcl(1) value of the dirichlet
!>                               - rcodcl(2) value of the exterior exchange
!>                                 coefficient (infinite if no exchange)
!>                               - rcodcl(3) value flux density
!>                                 (negative if gain) in w/m2 or roughness
!>                                 in m if icodcl=6
!>                                 -# for the velocity \f$ (\mu+\mu_T)
!>                                    \gradv \vect{u} \cdot \vect{n}  \f$
!>                                 -# for the pressure \f$ \Delta t
!>                                    \grad P \cdot \vect{n}  \f$
!>                                 -# for a scalar \f$ cp \left( K +
!>                                     \dfrac{K_T}{\sigma_T} \right)
!>                                     \grad T \cdot \vect{n} \f$
!> \param[in]     velipb        value of the velocity at \f$ \centip \f$
!>                               of boundary cells
!> \param[in]     rijipb        value of \f$ R_{ij} \f$ at \f$ \centip \f$
!>                               of boundary cells
!_______________________________________________________________________________


subroutine clsyvt &
 ( nscal  , icodcl ,                                     &
   rcodcl ,                                              &
   velipb , rijipb )

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

!===============================================================================
! Module files
!===============================================================================

use paramx
use numvar
use optcal
use cstphy
use cstnum
use pointe
use entsor
use albase
use field
use mesh
use cs_c_bindings

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

implicit none

! Arguments

integer          nscal

integer, pointer, dimension(:,:) :: icodcl

double precision, pointer, dimension(:,:,:) :: rcodcl
double precision, dimension(:,:) :: velipb
double precision, pointer, dimension(:,:) :: rijipb

! Local variables

integer          ifac, ii, isou
integer          iel, f_dim, clsyme
integer          iscal , ivar, idftnp
integer          kturt, turb_flux_model, turb_flux_model_type

double precision rnx, rny, rnz
double precision upx, upy, upz, usn
double precision tx, ty, tz, txn, t2x, t2y, t2z
double precision eloglo(3,3), alpha(6,6)
double precision srfbnf, rcodcn, hint, visclc, visctc, distbf
double precision cpp
double precision hintv(6), rnn(6), htnn(6)
double precision visci(3,3), fikis, viscis, distfi
double precision fcoefa(6), fcoefb(6), fcofaf(6), fcofbf(6), fcofad(6), fcofbd(6)

double precision, dimension(:,:), pointer :: coefau, cofafu, cfaale, claale
double precision, dimension(:,:), pointer :: visten
double precision, dimension(:,:,:), pointer :: coefbu, cofbfu, cfbale, clbale
double precision, dimension(:,:), pointer :: coefa_rij, coefaf_rij, coefad_rij
double precision, dimension(:,:,:), pointer :: coefb_rij, coefbf_rij, coefbd_rij

double precision, dimension(:), pointer :: viscl, visct
double precision, dimension(:), pointer :: cpro_visma_s
double precision, dimension(:,:), pointer :: cpro_visma_v

type(var_cal_opt) :: vcopt_uma, vcopt_rij

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

!===============================================================================
! 1. Initializations
!===============================================================================

cpp = 0.d0

if (itytur.eq.3 .and. idirsm.eq.1) call field_get_val_v(ivsten, visten)

call field_get_val_s(iviscl, viscl)
call field_get_val_s(ivisct, visct)

! Boundary Conditions

call field_get_coefa_v(ivarfl(iu), coefau)
call field_get_coefb_v(ivarfl(iu), coefbu)
call field_get_coefaf_v(ivarfl(iu), cofafu)
call field_get_coefbf_v(ivarfl(iu), cofbfu)

if (itytur.eq.3) then
  call field_get_key_struct_var_cal_opt(ivarfl(irij), vcopt_rij)
  call field_get_coefa_v(ivarfl(irij), coefa_rij)
  call field_get_coefb_v(ivarfl(irij), coefb_rij)
  call field_get_coefaf_v(ivarfl(irij), coefaf_rij)
  call field_get_coefbf_v(ivarfl(irij), coefbf_rij)
  call field_get_coefad_v(ivarfl(irij), coefad_rij)
  call field_get_coefbd_v(ivarfl(irij), coefbd_rij)
endif

! Get the turbulent flux model key id
call field_get_key_id('turbulent_flux_model', kturt)

! --- Begin the loop over boundary faces
do ifac = 1, nfabor

  ! --- Test sur la presence d'une condition de symetrie vitesse : debut
  if (icodcl(ifac,iu).eq.4) then

    ! --- To cancel the mass flux
    isympa(ifac) = 0

    ! Geometric quantities
    srfbnf = surfbn(ifac)

    !===========================================================================
    ! 1. Local framework
    !===========================================================================

    ! Unit normal

    rnx = surfbo(1,ifac)/srfbnf
    rny = surfbo(2,ifac)/srfbnf
    rnz = surfbo(3,ifac)/srfbnf

    ! En ALE, on a eventuellement une vitesse de deplacement de la face
    !   donc seule la composante normale importe (on continue a determiner
    !   TX a partir de la vitesse tangentielle absolue car l'orientation
    !   de TX et T2X est sans importance pour les symetries)
    rcodcn = 0.d0
    if (iale.ge.1) then
      rcodcn = rcodcl(ifac,iu,1)*rnx                           &
             + rcodcl(ifac,iv,1)*rny                           &
             + rcodcl(ifac,iw,1)*rnz
    endif

    upx = velipb(ifac,1)
    upy = velipb(ifac,2)
    upz = velipb(ifac,3)

    if (itytur.eq.3) then

      ! Relative tangential velocity

      usn = upx*rnx+upy*rny+upz*rnz
      tx  = upx -usn*rnx
      ty  = upy -usn*rny
      tz  = upz -usn*rnz
      txn = sqrt( tx**2 +ty**2 +tz**2 )

      ! Unit tangent

      if( txn.ge.epzero) then

        tx = tx/txn
        ty = ty/txn
        tz = tz/txn

      else

        ! If the velocity is zero,
        !  Tx, Ty, Tz is not used (we cancel the velocity), so we assign any
        !  value (zero for example)

        tx  = 0.d0
        ty  = 0.d0
        tz  = 0.d0

      endif

      ! --> T2 = RN X T (where X is the cross product)

      t2x = rny*tz - rnz*ty
      t2y = rnz*tx - rnx*tz
      t2z = rnx*ty - rny*tx

      ! --> Orthogonal matrix for change of reference frame ELOGLOij
      !     (from local to global reference frame)

      !                      |TX  -RNX  T2X|
      !             ELOGLO = |TY  -RNY  T2Y|
      !                      |TZ  -RNZ  T2Z|

      !    Its transpose ELOGLOt is its inverse

      eloglo(1,1) =  tx
      eloglo(1,2) = -rnx
      eloglo(1,3) =  t2x
      eloglo(2,1) =  ty
      eloglo(2,2) = -rny
      eloglo(2,3) =  t2y
      eloglo(3,1) =  tz
      eloglo(3,2) = -rnz
      eloglo(3,3) =  t2z

      ! Compute Reynolds stress transformation matrix

      clsyme = 1
      call turbulence_bc_rij_transform(clsyme, eloglo, alpha)

    endif

    !===========================================================================
    ! 2. Boundary conditions on the velocity
    !    (totaly implicit)
    !    The condition is a zero (except in ALE) Dirichlet on the normal component
    !                     a homogenous Neumann on the other components
    !===========================================================================

    ! Coupled solving of the velocity components

    iel = ifabor(ifac)
    ! --- Physical properties
    visclc = viscl(iel)
    visctc = visct(iel)

    ! --- Geometrical quantity
    distbf = distb(ifac)

    if (itytur.eq.3) then
      hint =   visclc         /distbf
    else
      hint = ( visclc+visctc )/distbf
    endif

    ! Gradient BCs
    coefau(1,ifac) = rcodcn*rnx
    coefau(2,ifac) = rcodcn*rny
    coefau(3,ifac) = rcodcn*rnz

    coefbu(1,1,ifac) = 1.d0-rnx**2
    coefbu(2,2,ifac) = 1.d0-rny**2
    coefbu(3,3,ifac) = 1.d0-rnz**2

    coefbu(1,2,ifac) = -rnx*rny
    coefbu(1,3,ifac) = -rnx*rnz
    coefbu(2,1,ifac) = -rny*rnx
    coefbu(2,3,ifac) = -rny*rnz
    coefbu(3,1,ifac) = -rnz*rnx
    coefbu(3,2,ifac) = -rnz*rny

    ! Flux BCs
    cofafu(1,ifac) = -hint*rcodcn*rnx
    cofafu(2,ifac) = -hint*rcodcn*rny
    cofafu(3,ifac) = -hint*rcodcn*rnz

    cofbfu(1,1,ifac) = hint*rnx**2
    cofbfu(2,2,ifac) = hint*rny**2
    cofbfu(3,3,ifac) = hint*rnz**2

    cofbfu(1,2,ifac) = hint*rnx*rny
    cofbfu(1,3,ifac) = hint*rnx*rnz
    cofbfu(2,1,ifac) = hint*rny*rnx
    cofbfu(2,3,ifac) = hint*rny*rnz
    cofbfu(3,1,ifac) = hint*rnz*rnx
    cofbfu(3,2,ifac) = hint*rnz*rny

    !===========================================================================
    ! 3. Boundary conditions on Rij (partially implicited)
    !===========================================================================

    if (itytur.eq.3) then

      ! Symmetric tensor diffusivity (Daly Harlow -- GGDH)
      if (iand(vcopt_rij%idften, ANISOTROPIC_RIGHT_DIFFUSION).ne.0) then

        visci(1,1) = visclc + visten(1,iel)
        visci(2,2) = visclc + visten(2,iel)
        visci(3,3) = visclc + visten(3,iel)
        visci(1,2) =          visten(4,iel)
        visci(2,1) =          visten(4,iel)
        visci(2,3) =          visten(5,iel)
        visci(3,2) =          visten(5,iel)
        visci(1,3) =          visten(6,iel)
        visci(3,1) =          visten(6,iel)

        ! ||Ki.S||^2
        viscis = ( visci(1,1)*surfbo(1,ifac)       &
                 + visci(1,2)*surfbo(2,ifac)       &
                 + visci(1,3)*surfbo(3,ifac))**2   &
               + ( visci(2,1)*surfbo(1,ifac)       &
                 + visci(2,2)*surfbo(2,ifac)       &
                 + visci(2,3)*surfbo(3,ifac))**2   &
               + ( visci(3,1)*surfbo(1,ifac)       &
                 + visci(3,2)*surfbo(2,ifac)       &
                 + visci(3,3)*surfbo(3,ifac))**2

        ! IF.Ki.S
        fikis = ( (cdgfbo(1,ifac)-xyzcen(1,iel))*visci(1,1)   &
                + (cdgfbo(2,ifac)-xyzcen(2,iel))*visci(2,1)   &
                + (cdgfbo(3,ifac)-xyzcen(3,iel))*visci(3,1)   &
                )*surfbo(1,ifac)                              &
              + ( (cdgfbo(1,ifac)-xyzcen(1,iel))*visci(1,2)   &
                + (cdgfbo(2,ifac)-xyzcen(2,iel))*visci(2,2)   &
                + (cdgfbo(3,ifac)-xyzcen(3,iel))*visci(3,2)   &
                )*surfbo(2,ifac)                              &
              + ( (cdgfbo(1,ifac)-xyzcen(1,iel))*visci(1,3)   &
                + (cdgfbo(2,ifac)-xyzcen(2,iel))*visci(2,3)   &
                + (cdgfbo(3,ifac)-xyzcen(3,iel))*visci(3,3)   &
                )*surfbo(3,ifac)

        distfi = distb(ifac)

        ! Take I" so that I"F= eps*||FI||*Ki.n when J" is in cell rji
        ! NB: eps =1.d-1 must be consistent with vitens.f90
        fikis = max(fikis, 1.d-1*sqrt(viscis)*distfi)

        hint = viscis/surfbn(ifac)/fikis

      ! Scalar diffusivity
      else
        hint = (visclc+visctc*csrij/cmu)/distbf
      endif

      do isou = 1, 6
        fcoefa(isou) = 0.d0
        fcofad(isou) = 0.d0
        fcoefb(isou) = 0.d0
        fcofbd(isou) = 0.d0
      enddo

      do isou = 1,6

        if (isou.eq.1) then
          ivar = ir11
        elseif (isou.eq.2) then
          ivar = ir22
        elseif (isou.eq.3) then
          ivar = ir33
        elseif (isou.eq.4) then
          ivar = ir12
        elseif (isou.eq.5) then
          ivar = ir23
        elseif (isou.eq.6) then
          ivar = ir13
        endif

        ! Partial (or total if coupled) implicitation
        if (irijco.eq.1) then
          coefa_rij(isou, ifac)  = 0.d0
          coefaf_rij(isou, ifac) = 0.d0
          coefad_rij(isou, ifac) = 0.d0
          do ii = 1, 6
            coefb_rij(isou,ii, ifac)  = alpha(ii, isou)
            if (ii.eq.isou) then
              coefbf_rij(isou,ii, ifac) = hint * (1.d0 - coefb_rij(isou,ii, ifac))
            else
              coefbf_rij(isou,ii, ifac) = - hint * coefb_rij(isou,ii, ifac)
            endif
            coefbd_rij(isou,ii, ifac) = coefb_rij(isou,ii, ifac)
          enddo

        else if (iclsyr.eq.1) then
          do ii = 1, 6
            if (ii.ne.isou) then
              fcoefa(isou) = fcoefa(isou) + alpha(isou,ii) * rijipb(ifac,ii)
            endif
          enddo
          fcoefb(isou) = alpha(isou,isou)
        else
          do ii = 1, 6
            fcoefa(isou) = fcoefa(isou) + alpha(isou,ii) * rijipb(ifac,ii)
          enddo
          fcoefb(isou) = 0.d0
        endif

        ! Boundary conditions for the momentum equation
        fcofad(isou) = fcoefa(isou)
        fcofbd(isou) = fcoefb(isou)

        ! Translate into Diffusive flux BCs
        fcofaf(isou) = -hint*fcoefa(isou)
        fcofbf(isou) = hint*(1.d0-fcoefb(isou))
      enddo

      if (irijco.ne.1) then
        do isou = 1, 6
          coefa_rij(isou,ifac) = fcoefa(isou)
          coefaf_rij(isou,ifac) = fcofaf(isou)
          coefad_rij(isou,ifac) = fcofad(isou)
          do ii = 1,6
            coefb_rij(isou,ii,ifac) = 0
            coefbf_rij(isou,ii,ifac) = 0
            coefbd_rij(isou,ii,ifac) = 0
          enddo
          coefb_rij(isou,isou,ifac) = fcoefb(isou)
          coefbf_rij(isou,isou,ifac) = fcofbf(isou)
          coefbd_rij(isou,isou,ifac) = fcofbd(isou)
        enddo
      endif

    endif

  endif
! --- Test sur la presence d'une condition de symetrie vitesse : fin

enddo
! ---  End of loop over boundary faces

!===========================================================================
! 3.bis Boundary conditions on transported vectors
!===========================================================================

do iscal = 1, nscal

  ! Get the associated turbulent flux model
  call field_get_key_int(ivarfl(isca(iscal)), kturt, turb_flux_model)
  turb_flux_model_type = turb_flux_model / 10

  ! u'T'
  if (turb_flux_model_type.eq.3) then
    call clsyvt_scalar(iscal, icodcl)
  endif

  ! additional transported vectors
  call field_get_dim(ivarfl(isca(iscal)), f_dim)
  if (f_dim.gt.1) then
    call clsyvt_vector(iscal, icodcl)
  endif
enddo

!===============================================================================
! 4. Symmetry boundary conditions for mesh velocity (ALE module)
!===============================================================================

if (iale.eq.1) then

  call field_get_coefa_v(ivarfl(iuma), claale)
  call field_get_coefb_v(ivarfl(iuma), clbale)
  call field_get_coefaf_v(ivarfl(iuma), cfaale)
  call field_get_coefbf_v(ivarfl(iuma), cfbale)

  call field_get_key_struct_var_cal_opt(ivarfl(iuma), vcopt_uma)
  idftnp = vcopt_uma%idften

  if (iand(idftnp, ISOTROPIC_DIFFUSION).ne.0) then
    call field_get_val_s(ivisma, cpro_visma_s)
  else if (iand(idftnp, ANISOTROPIC_LEFT_DIFFUSION).ne.0) then
    call field_get_val_v(ivisma, cpro_visma_v)
  endif

  do ifac = 1, nfabor
    if (icodcl(ifac,iuma).eq.4) then

      iel = ifabor(ifac)

      ! For a sliding boundary, the normal velocity is enforced to zero
      ! whereas the other components have an Homogenous Neumann
      ! NB: no recontruction in I' here

      ! --- Geometrical quantity
      distbf = distb(ifac)
      srfbnf = surfbn(ifac)

      ! --- Physical properties
      if (iand(idftnp, ISOTROPIC_DIFFUSION).ne.0) then
        do ii = 1, 3
          hintv(ii) = cpro_visma_s(iel)/distbf
          hintv(ii+3) = 0.d0
        enddo
      else if (iand(idftnp, ANISOTROPIC_LEFT_DIFFUSION).ne.0) then
        do ii = 1, 6
          hintv(ii) = cpro_visma_v(ii,iel)/distbf
        enddo
      endif

      ! Unit normal
      rnx = surfbo(1,ifac)/srfbnf
      rny = surfbo(2,ifac)/srfbnf
      rnz = surfbo(3,ifac)/srfbnf

      ! Coupled solving of the velocity components

      ! Gradient BCs
      claale(1,ifac) = 0.d0
      claale(2,ifac) = 0.d0
      claale(3,ifac) = 0.d0

      clbale(1,1,ifac) = 1.d0-rnx**2
      clbale(2,2,ifac) = 1.d0-rny**2
      clbale(3,3,ifac) = 1.d0-rnz**2

      clbale(1,2,ifac) = -rnx*rny
      clbale(2,1,ifac) = -rny*rnx
      clbale(1,3,ifac) = -rnx*rnz
      clbale(3,1,ifac) = -rnz*rnx
      clbale(2,3,ifac) = -rny*rnz
      clbale(3,2,ifac) = -rnz*rny

      ! Flux BCs
      cfaale(1,ifac) = 0.d0
      cfaale(2,ifac) = 0.d0
      cfaale(3,ifac) = 0.d0

      rnn(1) = rnx**2
      rnn(2) = rny**2
      rnn(3) = rnz**2
      rnn(4) = rnx*rny
      rnn(5) = rny*rnz
      rnn(6) = rnx*rnz

      call symmetric_matrix_product(hintv, rnn, htnn)
      cfbale(1,1,ifac) = htnn(1)
      cfbale(2,2,ifac) = htnn(2)
      cfbale(3,3,ifac) = htnn(3)
      cfbale(1,2,ifac) = htnn(4)
      cfbale(2,1,ifac) = htnn(4)
      cfbale(1,3,ifac) = htnn(6)
      cfbale(3,1,ifac) = htnn(6)
      cfbale(2,3,ifac) = htnn(5)
      cfbale(3,2,ifac) = htnn(5)

    endif
  enddo

endif

!----
! End
!----

return
end subroutine

!===============================================================================
! Local functions
!===============================================================================

!-------------------------------------------------------------------------------
! Arguments
!______________________________________________________________________________.
!  mode           name          role                                           !
!______________________________________________________________________________!
!> \param[in]     iscal         scalar id
!> \param[in,out] icodcl        face boundary condition code:
!>                               - 1 Dirichlet
!>                               - 3 Neumann
!>                               - 4 sliding and
!>                                 \f$ \vect{u} \cdot \vect{n} = 0 \f$
!>                               - 5 smooth wall and
!>                                 \f$ \vect{u} \cdot \vect{n} = 0 \f$
!>                               - 6 rough wall and
!>                                 \f$ \vect{u} \cdot \vect{n} = 0 \f$
!>                               - 9 free inlet/outlet
!>                                 (input mass flux blocked to 0)
!_______________________________________________________________________________

subroutine clsyvt_scalar &
 ( iscal  , icodcl )

!===============================================================================
! Module files
!===============================================================================

use paramx
use numvar
use optcal
use cstphy
use cstnum
use dimens, only: nvar
use pointe
use entsor
use albase
use parall
use ppppar
use ppthch
use ppincl
use cplsat
use mesh
use field

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

implicit none

! Arguments

integer          iscal

integer          icodcl(nfabor,nvar)

! Local variables

integer          ivar, f_id
integer          ifac, iel, isou, jsou
integer          iscacp, ifcvsl
integer          kturt, turb_flux_model, turb_flux_model_type

double precision cpp, rkl, visclc, visls_0
double precision distbf, srfbnf
double precision rnx, rny, rnz
double precision hintt(6)
double precision hint, qimp

character(len=80) :: fname

double precision, dimension(:), pointer :: val_s, crom
double precision, dimension(:), pointer :: coefap, coefbp, cofafp, cofbfp
double precision, dimension(:,:), pointer :: coefaut, cofafut, cofarut, visten
double precision, dimension(:,:,:), pointer :: coefbut, cofbfut, cofbrut
double precision, dimension(:), pointer :: a_al, b_al, af_al, bf_al

double precision, dimension(:), pointer :: viscl, viscls, cpro_cp

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

! Get the turbulent flux model for the scalar
call field_get_key_id('turbulent_flux_model', kturt)
call field_get_key_int(ivarfl(isca(iscal)), kturt, turb_flux_model)
turb_flux_model_type = turb_flux_model / 10

if (turb_flux_model_type.ne.3) return

ivar = isca(iscal)
f_id = ivarfl(ivar)

call field_get_val_s(ivarfl(ivar), val_s)
call field_get_val_v(ivsten, visten)
call field_get_val_s(iviscl, viscl)

call field_get_coefa_s(f_id, coefap)
call field_get_coefb_s(f_id, coefbp)
call field_get_coefaf_s(f_id, cofafp)
call field_get_coefbf_s(f_id, cofbfp)

call field_get_val_s(icrom, crom)

if (icp.ge.0) then
  call field_get_val_s(icp, cpro_cp)
endif

! Reference diffusivity of the primal scalar
call field_get_key_double(f_id, kvisl0, visls_0)

call field_get_key_int (f_id, kivisl, ifcvsl)
if (ifcvsl .ge. 0) then
  call field_get_val_s(ifcvsl, viscls)
endif

! Does the scalar behave as a temperature ?
call field_get_key_int(f_id, kscacp, iscacp)

! Turbulent diffusive flux of the scalar T
! (blending factor so that the component v'T' have only
!  mu_T/(mu+mu_T)* Phi_T)

! Name of the scalar ivar
call field_get_name(ivarfl(ivar), fname)

! Index of the corresponding turbulent flux
call field_get_id(trim(fname)//'_turbulent_flux', f_id)

call field_get_coefa_v(f_id,coefaut)
call field_get_coefb_v(f_id,coefbut)
call field_get_coefaf_v(f_id,cofafut)
call field_get_coefbf_v(f_id,cofbfut)
call field_get_coefad_v(f_id,cofarut)
call field_get_coefbd_v(f_id,cofbrut)

! EB-GGDH/AFM/DFM alpha boundary conditions
if (turb_flux_model.eq.11 .or. turb_flux_model.eq.21 .or. turb_flux_model.eq.31) then

  ! Name of the scalar ivar
  call field_get_name(ivarfl(ivar), fname)

  ! Index of the corresponding turbulent flux
  call field_get_id(trim(fname)//'_alpha', f_id)

  call field_get_coefa_s (f_id, a_al)
  call field_get_coefb_s (f_id, b_al)
  call field_get_coefaf_s(f_id, af_al)
  call field_get_coefbf_s(f_id, bf_al)
endif

! --- Loop on boundary faces
do ifac = 1, nfabor

  ! Test on symmetry boundary condition: start
  if (icodcl(ifac,iu).eq.4) then

    iel = ifabor(ifac)
    ! --- Physical Properties
    visclc = viscl(iel)
    cpp = 1.d0
    if (iscacp.eq.1) then
      if (icp.ge.0) then
        cpp = cpro_cp(iel)
      else
        cpp = cp0
      endif
    endif

    ! --- Geometrical quantities
    distbf = distb(ifac)
    srfbnf = surfbn(ifac)

    rnx = surfbo(1,ifac)/srfbnf
    rny = surfbo(2,ifac)/srfbnf
    rnz = surfbo(3,ifac)/srfbnf

    if (ifcvsl .lt. 0) then
      rkl = visls_0/cpp
    else
      rkl = viscls(iel)/cpp
    endif

    !FIXME for EB DFM
    do isou = 1, 6
      if (isou.le.3) then
        hintt(isou) = (0.5d0*(visclc + rkl)    &
          + ctheta(iscal)*visten(isou,iel)/csrij) / distbf
      else
        hintt(isou) = ctheta(iscal)*visten(isou,iel) / csrij / distbf
      endif
    enddo

    ! Gradient BCs
    coefaut(1,ifac) = 0.d0
    coefaut(2,ifac) = 0.d0
    coefaut(3,ifac) = 0.d0

    coefbut(1,1,ifac) = 1.d0-rnx**2
    coefbut(2,2,ifac) = 1.d0-rny**2
    coefbut(3,3,ifac) = 1.d0-rnz**2

    coefbut(1,2,ifac) = -rnx*rny
    coefbut(1,3,ifac) = -rnx*rnz
    coefbut(2,1,ifac) = -rny*rnx
    coefbut(2,3,ifac) = -rny*rnz
    coefbut(3,1,ifac) = -rnz*rnx
    coefbut(3,2,ifac) = -rnz*rny

    ! Flux BCs
    cofafut(1,ifac) = 0.d0
    cofafut(2,ifac) = 0.d0
    cofafut(3,ifac) = 0.d0

    cofbfut(1,1,ifac) = hintt(1)*rnx**2  + hintt(4)*rnx*rny + hintt(6)*rnx*rnz
    cofbfut(2,2,ifac) = hintt(4)*rnx*rny + hintt(2)*rny**2  + hintt(5)*rny*rnz
    cofbfut(3,3,ifac) = hintt(6)*rnx*rnz + hintt(5)*rny*rnz + hintt(3)*rnz**2

    cofbfut(1,2,ifac) = hintt(1)*rnx*rny + hintt(4)*rny**2  + hintt(6)*rny*rnz
    cofbfut(2,1,ifac) = hintt(4)*rnx**2  + hintt(2)*rny*rnx + hintt(5)*rnx*rnz
    cofbfut(1,3,ifac) = hintt(1)*rnx*rnz + hintt(4)*rny*rnz + hintt(6)*rnz**2
    cofbfut(3,1,ifac) = hintt(6)*rnx**2  + hintt(5)*rny*rnx + hintt(3)*rnz*rnx
    cofbfut(2,3,ifac) = hintt(4)*rnx*rnz + hintt(2)*rny*rnz + hintt(5)*rnz**2
    cofbfut(3,2,ifac) = hintt(6)*rnx*rny + hintt(5)*rny**2  + hintt(3)*rnz*rny

    ! Boundary conditions for thermal transport equation
    do isou = 1, 3
      cofarut(isou,ifac) = coefaut(isou,ifac)
      do jsou =1, 3
        cofbrut(isou,jsou,ifac) = coefbut(isou,jsou,ifac)
      enddo
    enddo

    ! EB-GGDH/AFM/DFM alpha boundary conditions
    if (turb_flux_model.eq.11 .or. turb_flux_model.eq.21 .or. turb_flux_model.eq.31) then

      ! Dirichlet Boundary Condition
      !-----------------------------

      qimp = 0.d0

      hint = 1.d0/distbf

      call set_neumann_scalar &
        !====================
      ( a_al(ifac), af_al(ifac),             &
        b_al(ifac), bf_al(ifac),             &
        qimp      , hint       )

    endif


  endif
  ! Test on velocity symmetry condition: end

enddo

!----
! End
!----

return
end subroutine

!===============================================================================
! Local functions
!===============================================================================

!-------------------------------------------------------------------------------
! Arguments
!______________________________________________________________________________.
!  mode           name          role                                           !
!______________________________________________________________________________!
!> \param[in]     iscal         scalar id
!> \param[in,out] icodcl        face boundary condition code:
!>                               - 1 Dirichlet
!>                               - 3 Neumann
!>                               - 4 sliding and
!>                                 \f$ \vect{u} \cdot \vect{n} = 0 \f$
!>                               - 5 smooth wall and
!>                                 \f$ \vect{u} \cdot \vect{n} = 0 \f$
!>                               - 6 rough wall and
!>                                 \f$ \vect{u} \cdot \vect{n} = 0 \f$
!>                               - 9 free inlet/outlet
!>                                 (input mass flux blocked to 0)
!_______________________________________________________________________________

subroutine clsyvt_vector &
 ( iscal  , icodcl )

!===============================================================================
! Module files
!===============================================================================

use paramx
use numvar
use optcal
use cstphy
use cstnum
use dimens, only: nvar
use pointe
use entsor
use albase
use parall
use ppppar
use ppthch
use ppincl
use cplsat
use mesh
use field
use cs_c_bindings

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

implicit none

! Arguments

integer          iscal

integer          icodcl(nfabor,nvar)

! Local variables

integer          ivar, f_id
integer          ifac, iel
integer          ifcvsl

double precision rkl, visclc, visls_0
double precision distbf, srfbnf
double precision rnx, rny, rnz, temp
double precision hintt(6)
double precision turb_schmidt

double precision, dimension(:), pointer :: crom
double precision, dimension(:,:), pointer :: coefav, cofafv
double precision, dimension(:,:,:), pointer :: coefbv, cofbfv
double precision, dimension(:,:), pointer :: visten

double precision, dimension(:), pointer :: viscl, visct, viscls

type(var_cal_opt) :: vcopt

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

ivar = isca(iscal)
f_id = ivarfl(ivar)

call field_get_key_struct_var_cal_opt(f_id, vcopt)

if (iand(vcopt%idften, ANISOTROPIC_DIFFUSION).ne.0) then
  if (iturb.ne.32) then
    call field_get_val_v(ivsten, visten)
  else ! EBRSM and (GGDH or AFM)
    call field_get_val_v(ivstes, visten)
  endif
endif

call field_get_val_s(iviscl, viscl)
call field_get_val_s(ivisct, visct)

call field_get_coefa_v(f_id, coefav)
call field_get_coefb_v(f_id, coefbv)
call field_get_coefaf_v(f_id, cofafv)
call field_get_coefbf_v(f_id, cofbfv)

call field_get_val_s(icrom, crom)

call field_get_key_int (f_id, kivisl, ifcvsl)
if (ifcvsl .ge. 0) then
  call field_get_val_s(ifcvsl, viscls)
endif

! retrieve turbulent Schmidt value for current scalar
call field_get_key_double(ivarfl(ivar), ksigmas, turb_schmidt)

! Reference diffusivity
call field_get_key_double(f_id, kvisl0, visls_0)

! --- Loop on boundary faces
do ifac = 1, nfabor

  ! Test on symmetry boundary condition: start
  if (icodcl(ifac,ivar).eq.4) then

    iel = ifabor(ifac)
    ! --- Physical Properties
    visclc = viscl(iel)

    ! --- Geometrical quantities
    distbf = distb(ifac)
    srfbnf = surfbn(ifac)

    rnx = surfbo(1,ifac)/srfbnf
    rny = surfbo(2,ifac)/srfbnf
    rnz = surfbo(3,ifac)/srfbnf

    if (ifcvsl .lt. 0) then
      rkl = visls_0
    else
      rkl = viscls(iel)
    endif

    ! Isotropic diffusivity
    if (iand(vcopt%idften, ISOTROPIC_DIFFUSION).ne.0) then
      hintt(1) = (vcopt%idifft*max(visct(iel),zero)/turb_schmidt + rkl)/distbf
      hintt(2) = hintt(1)
      hintt(3) = hintt(1)
      hintt(4) = 0.d0
      hintt(5) = 0.d0
      hintt(6) = 0.d0
    ! Symmetric tensor diffusivity
    elseif (iand(vcopt%idften, ANISOTROPIC_DIFFUSION).ne.0) then
      temp = vcopt%idifft*ctheta(iscal)/csrij
      hintt(1) = (temp*visten(1,iel) + rkl)/distbf
      hintt(2) = (temp*visten(2,iel) + rkl)/distbf
      hintt(3) = (temp*visten(3,iel) + rkl)/distbf
      hintt(4) =  temp*visten(4,iel)       /distbf
      hintt(5) =  temp*visten(5,iel)       /distbf
      hintt(6) =  temp*visten(6,iel)       /distbf
    endif

    ! Gradient BCs
    coefav(1,ifac) = 0.d0
    coefav(2,ifac) = 0.d0
    coefav(3,ifac) = 0.d0

    coefbv(1,1,ifac) = 1.d0-rnx**2
    coefbv(2,2,ifac) = 1.d0-rny**2
    coefbv(3,3,ifac) = 1.d0-rnz**2

    coefbv(1,2,ifac) = -rnx*rny
    coefbv(1,3,ifac) = -rnx*rnz
    coefbv(2,1,ifac) = -rny*rnx
    coefbv(2,3,ifac) = -rny*rnz
    coefbv(3,1,ifac) = -rnz*rnx
    coefbv(3,2,ifac) = -rnz*rny

    ! Flux BCs
    cofafv(1,ifac) = 0.d0
    cofafv(2,ifac) = 0.d0
    cofafv(3,ifac) = 0.d0

    cofbfv(1,1,ifac) = hintt(1)*rnx**2  + hintt(4)*rnx*rny + hintt(6)*rnx*rnz
    cofbfv(2,2,ifac) = hintt(4)*rnx*rny + hintt(2)*rny**2  + hintt(5)*rny*rnz
    cofbfv(3,3,ifac) = hintt(6)*rnx*rnz + hintt(5)*rny*rnz + hintt(3)*rnz**2

    cofbfv(1,2,ifac) = hintt(1)*rnx*rny + hintt(4)*rny**2  + hintt(6)*rny*rnz
    cofbfv(2,1,ifac) = hintt(1)*rnx*rny + hintt(4)*rny**2  + hintt(6)*rny*rnz
    cofbfv(1,3,ifac) = hintt(1)*rnx*rnz + hintt(4)*rny*rnz + hintt(6)*rnz**2
    cofbfv(3,1,ifac) = hintt(1)*rnx*rnz + hintt(4)*rny*rnz + hintt(6)*rnz**2
    cofbfv(2,3,ifac) = hintt(4)*rnx*rnz + hintt(2)*rny*rnz + hintt(5)*rnz**2
    cofbfv(3,2,ifac) = hintt(4)*rnx*rnz + hintt(2)*rny*rnz + hintt(5)*rnz**2

  endif
  ! Test on velocity symmetry condition: end

enddo

!----
! End
!----

return
end subroutine