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

!                      Code_Saturne version 3.0.0
!                      --------------------------
! This file is part of Code_Saturne, a general-purpose CFD tool.
!
! Copyright (C) 1998-2013 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 cs_user_boundary_conditions.f90
!>
!> \brief User subroutine which fills boundary conditions arrays
!> (\c icodcl, \c rcodcl) for unknown variables.
!>
!>
!> \section intro Introduction
!>
!> Here one defines boundary conditions on a per-face basis.
!>
!> Boundary faces may be selected using the \ref getfbr subroutine.
!>
!> \code getfbr(string, nelts, eltlst) \endcode
!>  - string is a user-supplied character string containing selection criteria;
!>  - nelts is set by the subroutine. It is an integer value corresponding to
!>    the number of boundary faces verifying the selection criteria;
!>  - lstelt is set by the subroutine. It is an integer array of size nelts
!>    containing the list of boundary faces verifying the selection criteria.
!>
!> string may contain:
!>  - references to colors (ex.: 1, 8, 26, ...)
!>  - references to groups (ex.: inlet, group1, ...)
!>  - geometric criteria (ex. x < 0.1, y >= 0.25, ...)
!>
!> These criteria may be combined using logical operators (\c and,\c or) and
!> parentheses.
!>
!> \par Example
!> \code 1 and (group2 or group3) and y < 1 \endcode
!> will select boundary faces
!> of color 1, belonging to groups 'group2' or 'group3' and with face center
!> coordinate y less than 1.
!>
!> Operators priority, from highest to lowest:
!>  '( )' > 'not' > 'and' > 'or' > 'xor'
!>
!> Similarly, interior faces and cells can be identified using the \ref getfac
!> and \ref getcel subroutines (respectively). Their syntax are identical to
!> \ref getfbr syntax.
!>
!> For a more thorough description of the criteria syntax, see the user guide.
!>
!>
!> \section bc_types Boundary condition types
!>
!> Boundary conditions may be assigned in two ways.
!>
!>
!> \subsection std_bcs For "standard" boundary conditions:
!>
!> One defines a code in the \c itypfb
!> array (of dimensions number of boundary faces).
!> This code will then be used by a non-user subroutine to assign the
!> following conditions.
!> The available codes are:
!>  - \c ientre: Inlet
!>  - \c isolib: Free outlet
!>  - \c isymet: Symmetry
!>  - \c iparoi: Wall (smooth)
!>  - \c iparug: Rough wall
!>
!> These integers are defined elsewhere (in paramx.f90 module).
!> Their value is greater than or equal to 1 and less than or  equal to
!> ntypmx (value fixed in paramx.h)
!>
!> In addition, some values must be defined:
!>  - Inlet (more precisely, inlet/outlet with prescribed flow, as the flow
!>           may be prescribed as an outflow):
!>    - Dirichlet conditions on variables other than pressure are mandatory
!>      if the flow is incoming, optional if the flow is outgoing (the code
!>      assigns zero flux if no Dirichlet is specified); thus,
!>      at face \c ifac, for the variable \c ivar: \c rcodcl(ifac, ivar, 1)
!>
!>
!>  - Smooth wall: (= impermeable solid, with smooth friction)
!>    - Velocity value for sliding wall if applicable:
!>                  - \c rcodcl(ifac, iu, 1) = fluid velocity in the x direction
!>                  - \c rcodcl(ifac, iv, 1) = fluid velocity in the y direction
!>                  - \c rcodcl(ifac, iw, 1) = fluid velocity in the z direction
!>    - Specific code and prescribed temperature value at wall if applicable:
!>                  - \c icodcl(ifac, ivar)    = 5
!>                  - \c rcodcl(ifac, ivar, 1) = prescribed temperature
!>    - Specific code and prescribed flux value at wall if applicable:
!>                  - \c icodcl(ifac, ivar)    = 3
!>                  - \c rcodcl(ifac, ivar, 3) = prescribed flux
!>    .
!>    Note that the default condition for scalars (other than k and epsilon)
!>    is homogeneous Neumann.
!>
!>
!>  - Rough wall: (= impermeable solid, with rough friction)
!>    - Velocity value for sliding wall if applicable:
!>                  - \c rcodcl(ifac, iu, 1) = fluid velocity in the x direction
!>                  - \c rcodcl(ifac, iv, 1) = fluid velocity in the y direction
!>                  - \c rcodcl(ifac, iw, 1) = fluid velocity in the z direction
!>    - Value of the dynamic roughness height to specify in
!>                  - \c rcodcl(ifac, iu, 3)
!>    - Value of the scalar roughness height (if required) to specify in
!>                  - \c rcodcl(ifac, iv, 3) (values for iw are not used)
!>    - Specific code and prescribed temperature value at wall if applicable:
!>                  - \c icodcl(ifac, ivar)    = 6
!>                  - \c rcodcl(ifac, ivar, 1) = prescribed temperature
!>    - Specific code and prescribed flux value at rough wall, if applicable:
!>                  - \c icodcl(ifac, ivar)    = 3
!>                  - \c rcodcl(ifac, ivar, 3) = prescribed flux
!>    .
!>    Note that the default condition for scalars (other than k and epsilon)
!>    is homogeneous Neumann.
!>
!>  - Symmetry (= slip wall):
!>    - Nothing to specify
!>
!>  - Free outlet (more precisely free inlet/outlet with prescribed pressure)
!>    - Nothing to prescribe for pressure and velocity. For scalars and
!>      turbulent values, a Dirichlet value may optionally be specified.
!>      The behavior is as follows:
!>          - pressure is always handled as a Dirichlet condition
!>          - if the mass flow is inflowing:
!>              one retains the velocity at infinity
!>              Dirichlet condition for scalars and turbulent values
!>               (or zero flux if the user has not specified a
!>                Dirichlet value)
!>          - if the mass flow is outflowing:
!>              one prescribes zero flux on the velocity, the scalars,
!>              and turbulent values
!>    .
!>    Note that the pressure will be reset to p0 on the first free outlet
!>    face found.
!>
!>
!-------------------------------------------------------------------------------
! Arguments
!______________________________________________________________________________.
!  mode           name          role                                           !
!______________________________________________________________________________!
!> \param[in]     nvar          total number of variables
!> \param[in]     nscal         total number of scalars
!> \param[out]    icodcl        boundary condition code:
!>                               - 1 Dirichlet
!>                               - 2 Radiative outlet
!>                               - 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 rought wall and
!>                                 \f$ \vect{u} \cdot \vect{n} = 0 \f$
!>                               - 9 free inlet/outlet
!>                                 (input mass flux blocked to 0)
!> \param[in]     itrifb        indirection for boundary faces ordering
!> \param[in,out] itypfb        boundary face types
!> \param[out]    izfppp        boundary face zone number
!> \param[in]     dt            time step (per cell)
!> \param[in]     rtp, rtpa     calculated variables at cell centers
!> \param[in]                    (at current and previous time steps)
!> \param[in]     propce        physical properties at cell centers
!> \param[in]     propfa        physical properties at interior face centers
!> \param[in]     propfb        physical properties at boundary face centers
!> \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 roughtness
!>                                 in m if icodcl=6
!>                                 -# for the velocity \f$ (\mu+\mu_T)
!>                                    \gradt \, \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$
!_______________________________________________________________________________

subroutine cs_user_boundary_conditions &
 ( nvar   , nscal  ,                                              &
   icodcl , itrifb , itypfb , izfppp ,                            &
   dt     , rtp    , rtpa   , propce , propfa , propfb ,          &
   rcodcl )

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

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

use paramx
use numvar
use optcal
use cstphy
use cstnum
use entsor
use parall
use period
use ihmpre
use ppppar
use ppthch
use coincl
use cpincl
use ppincl
use ppcpfu
use atincl
use atsoil
use ctincl
use elincl
use cs_fuel_incl
use mesh

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

implicit none

! Arguments

integer          nvar   , nscal

integer          icodcl(nfabor,nvarcl)
integer          itrifb(nfabor), itypfb(nfabor)
integer          izfppp(nfabor)

double precision dt(ncelet), rtp(ncelet,*), rtpa(ncelet,*)
double precision propce(ncelet,*)
double precision propfa(nfac,*), propfb(nfabor,*)
double precision rcodcl(nfabor,nvarcl,3)

! Local variables

! INSERT_VARIABLE_DEFINITIONS_HERE

integer, allocatable, dimension(:) :: lstelt

! Local variables
integer          ifac, iel, ii, ivar
integer          ilelt, nlelt

double precision uin, vin, win
!===============================================================================


!===============================================================================
! Initialization
!===============================================================================

allocate(lstelt(nfabor))  ! temporary array for boundary faces selection

! INSERT_MAIN_CODE_HERE

call getfbr("INLET", nlelt, lstelt)
do ilelt = 1, nlelt
   !
    ifac = lstelt(ilelt)
    iel = ifabor(ifac)
    itypfb(ifac)=ientre
   ! Define the inlet velocities here
    icodcl(ifac,iu)     = 1
    rcodcl(ifac, iu, 1) = 0.50d0
    icodcl(ifac,iu)     = 1
    rcodcl(ifac, iv, 1) = 0.00d0
    icodcl(ifac,iu)     = 1
    rcodcl(ifac, iw, 1) = 0.00d0
   ! Define the inlet temperature
   ! icodcl(ifac, itemp)    = 3
   ! rcodcl(ifac, itemp, 1) = 325.0d0
    icodcl(ifac, isca(1))    = 1
    rcodcl(ifac, isca(1), 1) = 325.0d0
enddo

call getfbr("OUTLET", nlelt, lstelt)
do ilelt = 1, nlelt
   !
    ifac = lstelt(ilelt)
    iel = ifabor(ifac)
    itypfb(ifac) = isolib 
    icodcl(ifac, ipr) = 3
    rcodcl(ifac, ipr, 3) = 0.d0
enddo

call getfbr("WAL", nlelt, lstelt)
!==========
!do ilelt = 1, nlelt
   !
   ifac = lstelt(ilelt)
   itypfb(ifac) = iparoi
   ! Define the wall heat flux
!   icodcl(ifac, ivar)    = 3
!   rcodcl(ifac, ivar, 3) = 1000.0d0
   icodcl(ifac, isca(3))    = 3
   rcodcl(ifac, isca(3), 3) = 1000.0d0
!enddo


deallocate(lstelt)  ! temporary array for boundary faces selection

return
end subroutine cs_user_boundary_conditions