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

!                      Code_Saturne version 4.0.2
!                      --------------------------
! This file is part of Code_Saturne, a general-purpose CFD tool.
!
! Copyright (C) 1998-2015 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 cs_user_boundary_conditions_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, lstelt) \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.
!>
!>
!> \subsection nonstd_bcs For "non-standard" conditions:
!>
!> Other than (inlet, free outlet, wall, symmetry), one defines
!>  - on one hand, for each face:
!>    - an admissible \c itypfb value (i.e. greater than or equal to 1 and
!>      less than or equal to \c ntypmx; see its value in paramx.h).
!>      The values predefined in paramx.h:
!>      \c ientre, \c isolib, \c isymet, \c iparoi, \c iparug are in this range,
!>      and it is preferable not to assign one of these integers to \c itypfb
!>      randomly or in an inconsiderate manner. To avoid this, one may use
!>      \c iindef if one wish to avoid checking values in paramx.h. \c iindef
!>      is an admissible value to which no predefined boundary condition
!>      is attached.
!>      Note that the \c itypfb array is reinitialized at each time step to
!>      the non-admissible value of 0. If one forgets to modify \c itypfb for
!>      a given face, the code will stop.
!>
!>  - and on the other hand, for each face and each variable:
!>    - a code
!>                      - \c icodcl(ifac, ivar)
!>    - three real values
!>                      - \c rcodcl(ifac, ivar, 1)
!>                      - \c rcodcl(ifac, ivar, 2)
!>                      - \c rcodcl(ifac, ivar, 3)
!>
!> \anchor icodcl \anchor rcodcl
!> The value of \c icodcl is taken from the following:
!>  - 1: Dirichlet      (usable for any variable)
!>  - 3: Neumann        (usable for any variable)
!>  - 4: Symmetry       (usable only for the velocity and components of
!>                       the Rij tensor)
!>  - 5: Smooth wall    (usable for any variable except for pressure)
!>  - 6: Rough wall     (usable for any variable except for pressure)
!>  - 9: Free outlet    (usable only for velocity)
!>  - 13: Dirichlet for the advection operator and
!>        Neumann for the diffusion operator
!>
!> The values of the 3 \c rcodcl components are:
!>  - \c rcodcl(ifac, ivar, 1):
!>     - Dirichlet for the variable          if \c icodcl(ifac, ivar) = 1 or 13
!>     - Wall value (sliding velocity, temp) if \c icodcl(ifac, ivar) = 5
!>     .
!>     The dimension of \c rcodcl(ifac, ivar, 1) is that of the
!>     resolved variable, for instance:
!>        - U (velocity in m/s),
!>        - T (temperature in degrees)
!>        - H (enthalpy in J/kg)
!>        - F (passive scalar in -)
!>  - \c rcodcl(ifac, ivar, 2):
!>       "exterior" exchange coefficient (between the prescribed value
!>                        and the value at the domain boundary)
!>                        rinfin = infinite by default
!>     - For velocities U,                in kg/(m2 s):
!>        \c rcodcl(ifac, ivar, 2) =          (viscl+visct) / d
!>     - For the pressure P,              in  s/m:
!>        \c rcodcl(ifac, ivar, 2) =                     dt / d
!>     - For temperatures T,              in Watt/(m2 degres):
!>        \c rcodcl(ifac, ivar, 2) = Cp*(viscls+visct/sigmas) / d
!>     - For enthalpies H,                in kg /(m2 s):
!>        \c rcodcl(ifac, ivar, 2) =    (viscls+visct/sigmas) / d
!>     - For other scalars F              in:
!>        \c rcodcl(ifac, ivar, 2) =    (viscls+visct/sigmas) / d
!>            (d has the dimension of a distance in m)
!>
!>  - \c rcodcl(ifac, ivar, 3) if \c icodcl(ifac, ivar) = 3 or 13:
!>      Flux density (< 0 if gain, n outwards-facing normal)
!>     - For velocities U,                in kg/(m s2) = J:
!>        \c rcodcl(ifac, ivar, 3) =         -(viscl+visct) * (grad U).n
!>     - For pressure P,                  in kg/(m2 s):
!>        \c rcodcl(ifac, ivar, 3) =                    -dt * (grad P).n
!>     - For temperatures T,              in Watt/m2:
!>        \c rcodcl(ifac, ivar, 3) = -Cp*(viscls+visct/sigmas) * (grad T).n
!>     - For enthalpies H,                in Watt/m2:
!>        \c rcodcl(ifac, ivar, 3) = -(viscls+visct/sigmas) * (grad H).n
!>     - For other scalars F              in:
!>        \c rcodcl(ifac, ivar, 3) = -(viscls+visct/sigmas) * (grad F).n
!>
!>  - \c rcodcl(ifac, ivar, 3) if \c icodcl(ifac, ivar) = 6:
!>      Roughness for the rough wall law
!>     - For velocities U, dynamic roughness
!>         \c rcodcl(ifac, iu, 3) = roughd
!>     - For other scalars, thermal roughness
!>         \c rcodcl(ifac, iv, 3) = rough
!>
!>
!> Note that if the user assigns a value to \c itypfb equal to \c ientre, \c isolib,
!> \c isymet, \c iparoi, or \c iparug and does not modify \c icodcl (zero value by
!>  default), \c itypfb will define the boundary condition type.
!>
!> To the contrary, if the user prescribes \c icodcl(ifac, ivar) (nonzero),
!> the values assigned to \c rcodcl will be used for the considered face
!> and variable (if \c rcodcl values are not set, the default values will
!> be used for the face and variable, so:
!>                          - \c rcodcl(ifac, ivar, 1) = 0.d0
!>                          - \c rcodcl(ifac, ivar, 2) = rinfin
!>                          - \c rcodcl(ifac, ivar, 3) = 0.d0)
!>
!> Especially, one may have for example:
!>  - set \c itypfb(ifac) = \c iparoi which prescribes default wall
!>    conditions for all variables at face ifac,
!>  - and define IN ADDITION for variable ivar on this face specific
!>    conditions by specifying \c icodcl(ifac, ivar) and the 3 \c rcodcl values.
!>
!> The user may also assign to \c itypfb a value not equal to \c ientre, \c isolib,
!> \c isymet, \c iparoi, \c iparug, \c iindef but greater than or equal to 1 and less
!> than or equal to ntypmx (see values in param.h) to distinguish groups
!> or colors in other subroutines which are specific to the case and in
!> which itypfb is accessible. In this case though it will be necessary
!> to prescribe boundary conditions by assigning values to icodcl and to
!> the 3 \c rcodcl fields (as the value of \c itypfb will not be predefined in
!> the code).
!>
!>
!> \subsection comp_bcs Boundary condition types for compressible flows
!>
!> For compressible flows, only predefined boundary conditions may
!> be assigned among: \c iparoi, \c isymet, \c iesicf, \c isspcf, \c isopcf, \c iephcf, \c ieqhcf
!>
!>  - \c iparoi : standard wall
!>  - \c isymet : standard symmetry
!>
!>  - \c iesicf, \c isspcf, \c isopcf, \c iephcf, \c ieqhcf : inlet/outlet
!>
!> For inlets/outlets, we can prescribe
!> a value for turbulence and passive scalars in \c rcodcl(.,.,1)
!> for the case in which the mass flux is incoming. If this is not
!> done, a zero flux condition is applied.
!>
!> - \c iesicf: prescribed inlet/outlet (for example supersonic inlet)
!>           the user prescribes the velocity and all thermodynamic variables
!> - \c isspcf: supersonic outlet
!>           the user does not prescribe anything
!> - \c isopcf: subsonic outlet with prescribed pressure
!>           the user presribes the pressure
!> - \c iephcf: mixed inlet with prescribed total pressure and enthalpy
!>           the user prescribes the total pressure and total enthalpy
!> - \c ieqhcf: subsonic inlet with prescribed mass and enthalpy flow
!>           to be implemented
!>
!>
!> \subsection cons_rul Consistency rules
!>
!> A few consistency rules between \c icodcl codes for variables with
!> non-standard boundary conditions:
!>
!>  - Codes for velocity components must be identical
!>  - Codes for Rij components must be identical
!>  - If code (velocity or Rij) = 4
!>    one must have code (velocity and Rij) = 4
!>  - If code (velocity or turbulence) = 5
!>    one must have code (velocity and turbulence) = 5
!>  - If code (velocity or turbulence) = 6
!>    one must have code (velocity and turbulence) = 6
!>  - If scalar code (except pressure or fluctuations) = 5
!>    one must have velocity code = 5
!>  - If scalar code (except pressure or fluctuations) = 6
!>    one must have velocity code = 6
!>
!>
!> \remarks
!>   - Caution: to prescribe a flux (nonzero) to Rij, the viscosity to take
!>              into account is viscl even if visct exists
!>              (visct=rho cmu k2/epsilon)
!>   - One have the ordering array for boundary faces from the previous time
!>       step (except for the fist one, where \c itrifb has not been set yet).
!>   - The array of boundary face types \c itypfb has been reset before
!>       entering the subroutine.
!>
!>
!> \subsubsection cs_user_bc_cell_id Cell values of some variables
!>
!> Cell value field ids
!>
!> - Density:                        \c iprpfl(irom)
!> - Dynamic molecular viscosity:    \c iprpfl(iviscl)
!> - Turbulent viscosity:            \c iprpfl(ivisct)
!> - Specific heat:                  \c iprpfl(icp)
!> - Diffusivity(lambda):            \c field_get_key_int(ivarfl(isca(iscal)), &
!>                                      kivisl, ...)
!>
!>
!> \subsubsection fac_id Faces identification
!>
!> - Density:                               \c field id \c ibrom
!> - Boundary mass flux (for convecting \c ivar):
!>     field id \c iflmab
!>     using \c field_get_key_int(ivarfl(ivar), kbmasf, iflmab)
!> - For other values: take as an approximation the value in the adjacent cell
!>                     i.e. as above with \c iel = ifabor(ifac).
!-------------------------------------------------------------------------------

!-------------------------------------------------------------------------------
! 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 rough wall and
!>                                 \f$ \vect{u} \cdot \vect{n} = 0 \f$
!>                               - 9 free inlet/outlet
!>                                 (input mass flux blocked to 0)
!>                               - 13 Dirichlet for the advection operator and
!>                                    Neumann for the diffusion operator
!> \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,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)
!>                                    \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     ,                                                       &
   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
use field
use turbomachinery
use iso_c_binding
use cs_c_bindings

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

implicit none

! Arguments

integer          nvar   , nscal

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

double precision dt(ncelet)
double precision, dimension(:,:), pointer :: cvara_vel
double precision, dimension(:), pointer :: imasfl, cvara_scalt, cvara_pr
double precision rcodcl(nfabor,nvarcl,3)
double precision,parameter::c = 0.5d0 

! Local variables

! INSERT_VARIABLE_DEFINITIONS_HERE

integer, allocatable, dimension(:) :: lstelt

! added allocatable array
integer, allocatable, dimension(:) :: lsteltm
integer, allocatable, dimension(:) :: icnt
double precision, allocatable, dimension(:) :: xm, ym, xmg, ymg, rxmg, rymg
double precision, allocatable, dimension(:) :: uloc, vloc, wloc, tloc
double precision, allocatable, dimension(:) :: xtabfluxcyl, xtabfluxcyl2, xtabsurfl

double precision, allocatable, dimension(:) :: ulocg, vlocg, wlocg, tlocg
double precision, allocatable, dimension(:) :: xtabfluxcylg, xtabfluxcylg2, xtabsurflg

double precision, allocatable, dimension(:) :: rulocg, rvlocg, rwlocg, rtlocg
double precision, allocatable, dimension(:) :: rxtabfluxcylg, rxtabfluxcylg2, rxtabsurflg
! Local variables
integer          ifac, iel, ii, ivar, jj
integer          ilelt, nlelt, ileltm

double precision uin, vin, win, xin, yin, zin, courant

! added integer variables
integer         ielm1, ielm2, iscal
integer         nleltm, nleltg, nleltmg, inleltmg
integer         idx0, idx
integer         nlelt0, nleltm0
integer         iflmas
integer         i1, i2
! added double precision variables
double precision rxi, ryi, rxj, ryj
double precision xeps
double precision surf, ipond
double precision surfin, mmflux, mmflux2
double precision ubulk
!added characters
Character(10) :: nip
!===============================================================================


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

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

allocate(lsteltm(nfac))  ! temporary array for boundary faces selection
! INSERT_ADDITIONAL_INITIALIZATION_CODE_HERE
! Get the coordinates of the face centers at the mapped face.
call getfbr("INLET1", nlelt, lstelt)

iscal = iscalt         ! temperature scalar number
ivar =  isca(iscal)    ! temperature variable number

call field_get_val_prev_s(ivarfl(ivar), cvara_scalt)
call field_get_val_prev_s(ivarfl(ipr), cvara_pr)
call field_get_val_prev_v(ivarfl(iu), cvara_vel)

if (.true.) then! test whether the recycling need to be turned on.
  call getfac('INLET1A', nleltm, lsteltm)  ! select internal faces of 'INLET2A'.
                                         ! the returned number of elements 'nleltm'.
                                         ! the returned list of elements 'lsteltm'.

  deallocate(lsteltm)  ! temporary array for mapped faces selection

  nleltm0=nleltm
  nlelt0 =nlelt
  allocate(lsteltm(nleltm0))  ! temporary array for boundary faces selection

  call getfac('INLET1A', nleltm, lsteltm)  ! select internal faces of 'INLET2A'.

  ! (parallel action)global sum
  if (irangp.ge.0) then
    call parcpt(nleltm0)
    call parcpt(nlelt0)
  endif

!~   if (irangp.ge.0) then
!~       Write( nip, '(i10)' ) irangp
!~       open (unit = irangp, file = "lvar"//nip)
!~   endif
!~ 
!~ !  if (ntcabs .le. 3202) then
!~   if (irangp.ge.0) then
!~       write(irangp,*) 'nleltmg= ', nleltm0, 'irangp= ', irangp ! total number of elements on 'INLET1A'
!~       write(irangp,*) 'nleltg=', nlelt0, 'irangp= ', irangp    ! total number of elements on 'INLET1'
!~       write(irangp,*) 'nleltm=', nleltm, 'irangp= ', irangp    ! number of element on 'INLET1A' on local node
!~       write(irangp,*) 'nlelt=', nlelt, 'irangp= ', irangp      ! number of element on 'INLET1' on local node
!~       write(irangp,*) 'nfabor=', nfabor, 'irangp= ', irangp    ! number of boundary faces on the node
!~       write(irangp,*) 'nfab=', nfac, 'irangp= ', irangp        ! number of interior faces on the node
!~   endif
!  endif

  nleltmg = nleltm0
  nleltg = nlelt0

  allocate(xm(nleltm))                    ! temporary array for x-cor of the face centres of mapped faces selection on the local node.
  allocate(ym(nleltm))                    ! temporary array for y-cor of the face centres of mapped faces selection on the local node.

! Transfer the x and y coordinate to the user defined array xm() and ym(). The length of the arraies is 'nlelm'
  do ileltm = 1, nleltm
      ifac = lsteltm(ileltm)             ! index of internal faces.
      xm(ileltm) = cdgfac(1,ifac)        ! x-cor of the face centres of mapped faces selection.
      ym(ileltm) = cdgfac(2,ifac)        ! y-cor of the face centres of mapped faces selection.
  enddo

! (parallel action) send coordinate to global arrays
  allocate(xmg(nleltmg))	! temporary array for x-cor of the face centres of mapped faces selection on the global node.
  allocate(ymg(nleltmg))	! temporary array for y-cor of the face centres of mapped faces selection on the global node.
  if (irangp.ge.0) then
    call paragv(nleltm,nleltmg, xm, xmg)
    call paragv(nleltm,nleltmg, ym, ymg)
  endif
!  if (ntcabs .le. 3202) then
!~   if (irangp.ge.0) then
!~      write(irangp,*) 'maxxmg= ', maxval(xmg), 'maxymg= ', maxval(ymg)
!~      write(irangp,*) 'minxmg= ', minval(xmg), 'minymg= ', minval(ymg)
!~      write(irangp,*) 'maxxm= ', maxval(xm), 'maxym= ', maxval(ym)
!~      write(irangp,*) 'minxm= ', minval(xm), 'minym= ', minval(ym)
!~   endif
!  endif
! Get the velocity component from the mapped faces on the local nodes
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  allocate(uloc(nleltm))
  allocate(vloc(nleltm))
  allocate(wloc(nleltm))

  allocate(tloc(nleltm))

  allocate(xtabfluxcyl(nleltm))
  allocate(xtabfluxcyl2(nleltm))
  allocate(xtabsurfl(nleltm))
  
  call field_get_key_int(ivarfl(iw), kimasf, iflmas)
  call field_get_val_s(iflmas, imasfl)
!~   iflmas = ipprof(ifluma(iw)) !Rank i in propfa(.,i) of the properties defined at the internal faces.
!  if (ntcabs .le. 3202) then
!~   if (irangp.ge.0)then
!~      write(irangp,*) 'iflmas=', iflmas
!~   endif
!  endif

  do ileltm = 1, nleltm
     ifac = lsteltm(ileltm)
     i1 = ifacel(1,ifac)       ! Index-numbers of the 1st neighboring cells for each internal face.
     i2 = ifacel(2,ifac)       ! Index-numbers of the 1st neighboring cells for each internal face.
     ipond = pond(ifac)        ! pond(nfac), real array giving
                               !F J.n
                               !----- for every internal face. With regard to the mesh quality, its ideal value is 0.5.
                               !IJ.n
!   write(nfecra,*) 'ipond=', ipond   ! Test z-coordinate if the inlet boundary face.
     uloc(ileltm) = ipond*cvara_vel(1,i1)+(1.d0-ipond)*cvara_vel(1,i2) ! cvara_vel (ncelet,nvar) Array storing the values of the solved variables 
                                                               ! at the previous time step.
     vloc(ileltm) = ipond*cvara_vel(2,i1)+(1.d0-ipond)*cvara_vel(2,i2)
     wloc(ileltm) = ipond*cvara_vel(3,i1)+(1.d0-ipond)*cvara_vel(3,i2)

     tloc(ileltm) = ipond*cvara_scalt(i1)+(1.d0-ipond)*cvara_scalt(i2)

     xtabfluxcyl(ileltm) = imasfl(ifac)*surfac(3,ifac)/abs(surfac(3,ifac)) ! mass flow rate crossing the each selected internal face 
                                                                                  ! with direction normal to the faces.  
     surf = sqrt(surfac(1,ifac)**2+surfac(2,ifac)**2+surfac(3,ifac)**2)    ! area of the internal surface.
     xtabfluxcyl2(ileltm) = ro0*wloc(ileltm)*surf                          ! mass flow rate calculated by formula rho*w*A.
     xtabsurfl(ileltm) = surf		! array of area of each elements on the internal surface on local nodes.
  enddo

  ! (parallel action) send velocity, mass flow rate and aera of each faces to global arrays
  allocate(ulocg(nleltmg))		! allocate memory to the global array of x-velocity.
  allocate(vlocg(nleltmg))		! allocate memory to the global array of y-velocity.
  allocate(wlocg(nleltmg))		! allocate memory to the global array of z-velocity.

  allocate(tlocg(nleltmg))		! allocate memory to the global array of temperature.

  allocate(xtabfluxcylg(nleltmg))       ! allocate memory to the global array of mass flow rate crossing the each selected internal face.
  allocate(xtabfluxcylg2(nleltmg))      ! allocate memory to the global array of mass flow rate crossing the each selected internal face,     
					! based on fomula.
  allocate(xtabsurflg(nleltmg))		! allocate memory to the global array of area of each selected internal face.

  if (irangp.ge.0) then
    call paragv(nleltm,nleltmg, uloc, ulocg)
    call paragv(nleltm,nleltmg, vloc, vlocg)
    call paragv(nleltm,nleltmg, wloc, wlocg)

    call paragv(nleltm,nleltmg, tloc, tlocg)

    call paragv(nleltm,nleltmg, xtabfluxcyl, xtabfluxcylg)
    call paragv(nleltm,nleltmg, xtabfluxcyl2, xtabfluxcylg2)
    call paragv(nleltm,nleltmg, xtabsurfl, xtabsurflg)
  endif

  xeps = 1.0d-12
! (parallel action) mark the dupicated faces
  do ii = 1, nleltmg-1
     rxi=xmg(ii)
     ryi=ymg(ii)
     do jj = ii+1, nleltmg
        rxj = xmg(jj)
        ryj = ymg(jj)
        if (abs(rxi-rxj) .lt. xeps .and. abs(ryi-ryj) .lt. xeps) then
           xmg(jj)=1.0d12
           ymg(jj)=1.0d12
        endif
     enddo
  enddo

! (parallel action) find the number of the faces (avoid the deplicated faces)
  inleltmg = 0
  do ii = 1, nleltmg
      if (xmg(ii) .lt. 1.0d10 .and. ymg(ii) .lt. 1.0d10)then
         inleltmg = inleltmg+1
      endif
  enddo

!  if (ntcabs .le. 3202) then
!~   if (irangp.ge.0)then
!~      write(irangp,*) 'inleltmg=', inleltmg
!~   endif
!  endif

! (parallel action) set the un-duplicated coordinates information to 'rxmg' and 'rymg'
  allocate(rxmg(inleltmg))                    ! temporary array for x-cor of the face centres of mapped faces selection.
  allocate(rymg(inleltmg))                    ! temporary array for y-cor of the face centres of mapped faces selection.

  allocate(rulocg(inleltmg))                  ! temporary array for x-vel of the face centres of mapped faces selection.
  allocate(rvlocg(inleltmg))                  ! temporary array for y-vel of the face centres of mapped faces selection.
  allocate(rwlocg(inleltmg))                  ! temporary array for z-vel of the face centres of mapped faces selection.

  allocate(rtlocg(inleltmg))                  ! temporary array for temp of the face centres of mapped faces selection.

  allocate(rxtabfluxcylg(inleltmg))           ! temporary array for mass flow rate through the selected faces on mapped faces selection.
  allocate(rxtabfluxcylg2(inleltmg))          ! temporary array for mass flow rate through the selected faces on mapped faces selection.
  allocate(rxtabsurflg(inleltmg))             ! temporary array for area of select faces on mapped faces selection.

  idx0 = 0
  do ii = 1, nleltmg
      if (xmg(ii) .lt. 1.0d10 .and. ymg(ii) .lt. 1.0d10)then
         idx0 = idx0+1
         rxmg(idx0)=xmg(ii)
         rymg(idx0)=ymg(ii)

         rulocg(idx0)=ulocg(ii)
         rvlocg(idx0)=vlocg(ii)
         rwlocg(idx0)=wlocg(ii)
  
         rtlocg(idx0)=tlocg(ii)
  
         rxtabfluxcylg(idx0)=xtabfluxcylg(ii)
         rxtabfluxcylg2(idx0)=xtabfluxcylg2(ii)
         rxtabsurflg(idx0)=xtabsurflg(ii)
      endif
  enddo

!  if (ntcabs .le. 3202) then
!~   if (irangp.ge.0)then
!~      write(irangp,*) 'idx0=', idx0
!~   endif
!  endif

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Find the index number of mapped faces.
  allocate(icnt(nlelt)) ! The size of array icnt == number of inlet faces on the node
  idx = 0
  do ilelt = 1, nlelt
     ifac = lstelt(ilelt)
     xin = cdgfbo (1,ifac)
     yin = cdgfbo (2,ifac)
     zin = cdgfbo (3,ifac)
     do ileltm = 1, inleltmg
        if(abs(xin-rxmg(ileltm)) .le. xeps .and. abs(yin-rymg(ileltm)) .le. xeps) then
          idx = idx + 1
          icnt(idx) = ileltm   ! Transfer the index number of the mapped face to the array 'icnt'
        endif
     enddo
  enddo

!  if (ntcabs .le. 3202) then
!~   if (irangp.ge.0)then
!~      write(irangp,*) 'idx=', idx
!~   endif
!  endif
!!!!!!!!!!!!!!!!--------------------Upwards are for parallel02--------------------!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Calculate the global total mass flow rate from mapped face.
  mmflux =0.0d0
  mmflux2=0.0d0
  surfin =0.0d0
  do ii = 1, idx0
     mmflux  = mmflux + rxtabfluxcylg(ii)
     mmflux2 = mmflux2 + rxtabfluxcylg2(ii)
     surfin = surfin + rxtabsurflg(ii)
  enddo

!  if (irangp.ge.0) then
!    call parsom(mmflux)
!    call parsom(mmflux2)
!    call parsom(surfin)
!  endif

!~   if (irangp.ge.0)then
!~      write(irangp,*) 'mmflux=', mmflux
!~      write(irangp,*) 'mmflux2=', mmflux2
!~      write(irangp,*) 'surfin=', surfin
!~   endif

endif ! test of whether need to turning on recycling ended.

! INSERT_MAIN_CODE_HERE
!call getfbr("INLET1", nlelt, lstelt)
ubulk = -3.88d-2
!surfin = 0.00384d0

!~ if (irangp.ge.0)then
!~ !!   Write( nip, '(i10)' ) irangp                 !Just for uniform velocity distribution.
!~ !!   open (unit = irangp, file = "lvar"//nip)     !Just for uniform velocity distribution.
!~    write(irangp,*) 'mfluxin=', ubulk*ro0*surfin
!~ endif

if (ntcabs .le. 4000) then
  do ilelt = 1, nlelt
      ifac = lstelt(ilelt)
      iel = ifabor(ifac)
      itypfb(ifac)=ientre
      icodcl(ifac,iu)     = 1
      rcodcl(ifac, iu, 1) = 0.00d0
      icodcl(ifac,iv)     = 1
      rcodcl(ifac, iv, 1) = 0.00d0
      icodcl(ifac,iw)     = 1
      rcodcl(ifac, iw, 1) = ubulk*1.0d0
      icodcl(ifac, isca(1))    = 1
      rcodcl(ifac, isca(1), 1) = 16.97d5
  enddo
!~ elseif (ntcabs .gt. 4010 .and. ntcabs .le. 4050) then
!~     do ilelt = 1, nlelt
!~       ifac = lstelt(ilelt)
!~       iel = ifabor(ifac)
!~       itypfb(ifac)=ientre
!~       icodcl(ifac,iu)     = 1
!~       rcodcl(ifac, iu, 1) = rulocg(icnt(ilelt))*ubulk*5d0*ro0*surfin/mmflux
!~       icodcl(ifac,iu)     = 1
!~       rcodcl(ifac, iv, 1) = rvlocg(icnt(ilelt))*ubulk*5d0*ro0*surfin/mmflux
!~       icodcl(ifac,iu)     = 1
!~       rcodcl(ifac, iw, 1) = rwlocg(icnt(ilelt))*ubulk*5d0*ro0*surfin/mmflux
!~       icodcl(ifac, isca(1))    = 1
!~       rcodcl(ifac, isca(1), 1) = rtlocg(icnt(ilelt))
!~     enddo
elseif (ntcabs .gt. 4000 .and. ntcabs .le. 19200000) then
    do ilelt = 1, nlelt
      ifac = lstelt(ilelt)
      iel = ifabor(ifac)
      itypfb(ifac)=ientre
      icodcl(ifac,iu)     = 1
      rcodcl(ifac, iu, 1) = rulocg(icnt(ilelt))*ubulk*1d0*ro0*surfin/mmflux
      icodcl(ifac,iv)     = 1
      rcodcl(ifac, iv, 1) = rvlocg(icnt(ilelt))*ubulk*1d0*ro0*surfin/mmflux
      icodcl(ifac,iw)     = 1
      rcodcl(ifac, iw, 1) = rwlocg(icnt(ilelt))*ubulk*1d0*ro0*surfin/mmflux
      icodcl(ifac, isca(1))    = 1
      rcodcl(ifac, isca(1), 1) = rtlocg(icnt(ilelt))
     enddo
!~ elseif (ntcabs .gt. 19200 .and. ntcabs .le. 28800) then
!~     do ilelt = 1, nlelt
!~       ifac = lstelt(ilelt)
!~       iel = ifabor(ifac)
!~       itypfb(ifac)=ientre
!~       icodcl(ifac,iu)     = 1
!~       rcodcl(ifac, iu, 1) = rulocg(icnt(ilelt))*ubulk*1.0d0*ro0*surfin/mmflux
!~       icodcl(ifac,iu)     = 1
!~       rcodcl(ifac, iv, 1) = rvlocg(icnt(ilelt))*ubulk*1.0d0*ro0*surfin/mmflux
!~       icodcl(ifac,iu)     = 1
!~       rcodcl(ifac, iw, 1) = rwlocg(icnt(ilelt))*ubulk*1.0d0*ro0*surfin/mmflux
!~       icodcl(ifac, isca(1))    = 1
!~       rcodcl(ifac, isca(1), 1) = rtlocg(icnt(ilelt))
!~      enddo
endif

call getfbr("OUTLET4", nlelt, lstelt)
do ilelt = 1, nlelt
   !
   ifac = lstelt(ilelt)
   iel = ifabor(ifac)
   itypfb(ifac) = ifrent
   iel = ifabor(ifac)
   
   courant = (c)*dt(iel)/abs(distb(ifac))

   icodcl(ifac,ipr) = 2
   rcodcl(ifac,ipr,2) = courant
   rcodcl(ifac,ipr,1) = cvara_pr(iel)
   
   icodcl(ifac,iu) = 2
   rcodcl(ifac,iu,2) = courant
   rcodcl(ifac,iu,1) = cvara_vel(1,iel)

   icodcl(ifac,iv) = 2
   rcodcl(ifac,iv,2) = courant
   rcodcl(ifac,iv,1) = cvara_vel(2,iel)

   icodcl(ifac,iw) = 2
   rcodcl(ifac,iw,2) = courant
   rcodcl(ifac,iw,1) = cvara_vel(3,iel)
enddo

call getfbr("INNER-WALL1", nlelt, lstelt)
!==========
do ilelt = 1, nlelt
   !
   ifac = lstelt(ilelt)
   itypfb(ifac) = iparoi
   icodcl(ifac, isca(1))    = 3
   rcodcl(ifac, isca(1), 3) = 0.d0
enddo

call getfbr("INNER-WALL2", nlelt, lstelt)
!==========
do ilelt = 1, nlelt
   !
   ifac = lstelt(ilelt)
   itypfb(ifac) = iparoi
   icodcl(ifac, isca(1))    = 3
   rcodcl(ifac, isca(1), 3) = -1.0d4
enddo

call getfbr("INNER-WALL3", nlelt, lstelt)
!==========
do ilelt = 1, nlelt
   !
   ifac = lstelt(ilelt)
   itypfb(ifac) = iparoi
   icodcl(ifac, isca(1))    = 3
   rcodcl(ifac, isca(1), 3) = -1.0d4
enddo

call getfbr("INNER-WALL4", nlelt, lstelt)
!==========
do ilelt = 1, nlelt
   !
   ifac = lstelt(ilelt)
   itypfb(ifac) = iparoi
   icodcl(ifac, isca(1))    = 3
   rcodcl(ifac, isca(1), 3) = -1.0d4
enddo

call getfbr("OUTER-WALL1", nlelt, lstelt)
!==========
do ilelt = 1, nlelt
   !
   ifac = lstelt(ilelt)
   itypfb(ifac) = iparoi
   icodcl(ifac, isca(1))    = 3
   rcodcl(ifac, isca(1), 3) = 0.d0
enddo

call getfbr("OUTER-WALL2", nlelt, lstelt)
!==========
do ilelt = 1, nlelt
   !
   ifac = lstelt(ilelt)
   itypfb(ifac) = iparoi
   icodcl(ifac, isca(1))    = 3
   rcodcl(ifac, isca(1), 3) = 0.d0
enddo

call getfbr("OUTER-WALL3", nlelt, lstelt)
!==========
do ilelt = 1, nlelt
   !
   ifac = lstelt(ilelt)
   itypfb(ifac) = iparoi
   icodcl(ifac, isca(1))    = 3
   rcodcl(ifac, isca(1), 3) = 0.d0
enddo

call getfbr("OUTER-WALL4", nlelt, lstelt)
!==========
do ilelt = 1, nlelt
   !
   ifac = lstelt(ilelt)
   itypfb(ifac) = iparoi
   icodcl(ifac, isca(1))    = 3
   rcodcl(ifac, isca(1), 3) = 0.d0
enddo

!--------
! Formats
!--------

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

deallocate(xm) ! temporary array for x-cor of the face centres of mapped faces selection
deallocate(ym) ! temporary array for y-cor of the face centres of mapped faces selection
deallocate(xmg) ! temporary array for x-cor of the face centres of mapped faces selection
deallocate(ymg) ! temporary array for y-cor of the face centres of mapped faces selection
deallocate(uloc) ! temporary array for x-cor of the face centres of mapped faces selection
deallocate(vloc) ! temporary array for y-cor of the face centres of mapped faces selection
deallocate(wloc) ! temporary array for x-cor of the face centres of mapped faces selection
deallocate(xtabfluxcyl) ! temporary array for y-cor of the face centres of mapped faces selection
deallocate(xtabfluxcyl2) ! temporary array for y-cor of the face centres of mapped faces selection
deallocate(xtabsurfl) ! temporary array for y-cor of the face centres of mapped faces selection
deallocate(lstelt)  ! temporary array for boundary faces selection
deallocate(lsteltm)  ! temporary array for mapped faces selection
deallocate(tloc) ! temporary array for x-cor of the face centres of mapped faces selection
deallocate(ulocg)
deallocate(vlocg)
deallocate(wlocg)
deallocate(tlocg)
deallocate(xtabfluxcylg)
deallocate(xtabfluxcylg2)
deallocate(xtabsurflg)
deallocate(rxmg)                    ! temporary array for x-cor of the face centres of mapped faces selection.
deallocate(rymg)                    ! temporary array for y-cor of the face centres of mapped faces selection.
deallocate(rulocg)                    ! temporary array for x-cor of the face centres of mapped faces selection.
deallocate(rvlocg)                    ! temporary array for y-cor of the face centres of mapped faces selection.
deallocate(rwlocg)                    ! temporary array for y-cor of the face centres of mapped faces selection.
deallocate(rtlocg)                    ! temporary array for y-cor of the face centres of mapped faces selection.
deallocate(rxtabfluxcylg)                    ! temporary array for x-cor of the face centres of mapped faces selection.
deallocate(rxtabfluxcylg2)                    ! temporary array for y-cor of the face centres of mapped faces selection.
deallocate(rxtabsurflg)                    ! temporary array for y-cor of the face centres of mapped faces selection.

return
end subroutine cs_user_boundary_conditions

!> \section examples Examples
!>   Several examples are provided
!>   \ref cs_user_boundary_conditions_examples "here".