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

!                      Code_Saturne version 3.2.2
!                      --------------------------
! 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.

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




subroutine uslapr &
!================

 ( idvar  , iepart , izone  , iclass ,                            &
   nvar   , nscal  ,                                              &
   nbpmax , nvp    , nvp1   , nvep   , nivep  ,                   &
   ntersl , nvlsta , nvisbr ,                                     &
   itypfb , itrifb , itepa  , ifrlag ,                            &
   xxpart , yypart , zzpart ,                                     &
   tvpart , uupart , vvpart , wwpart , ddpart , ttpart  ,         &
   dt     , rtpa   , propce ,                                     &
   ettp   , tepa   )

!===============================================================================
! Purpose:
! ----------

!   Subroutine of the Lagrangian particle-tracking module:
!   -------------------------------------

!   User subroutine for the boundary conditions associated to
!   the particles (inlet and treatment of the other boundaries)
!
!   It allows to impose the values of the velocity, the diameter
!   and the temperature for the treated particle.
!
!   if idvar = 0 ==> the volume fraction is retrieved
!   if idvar = 1 ==> the 3 components of the velocity are retrieved
!   if idvar = 2 ==> the diameter is retrieved
!   if idvar = 3 ==> the temperature is retrieved


!-------------------------------------------------------------------------------
! Arguments
!__________________.____._____.________________________________________________.
! name             !type!mode ! role                                           !
!__________________!____!_____!________________________________________________!
! idvar            ! i  ! <-- ! type of the value(s) ta calculate              !
! iepart           ! i  ! <-- ! number of the particle cell                    !
! izone            ! i  ! <-- ! number of the particle zone                    !
! iclass           ! i  ! <-- ! number of the particle class                   !
! nvar             ! i  ! <-- ! total number of variables                      !
! nscal            ! i  ! <-- ! total number of scalars                        !
! nbpmax           ! i  ! <-- ! maximum number of particles allowed            !
! nvp              ! i  ! <-- ! number of particle variables                   !
! nvp1             ! i  ! <-- ! nvp minus position, fluid and part. velocities !
! nvep             ! i  ! <-- ! number of particle properties (integer)        !
! nivep            ! i  ! <-- ! number of particle properties (integer)        !
! ntersl           ! i  ! <-- ! number of source terms of return coupling      !
! nvlsta           ! i  ! <-- ! nb of Lagrangian statistical variables         !
! nvisbr           ! i  ! <-- ! number of boundary statistics                  !
! itrifb(nfabor)   ! ia ! <-- ! indirection for the sorting of the             !
! itypfb(nfabor)   ! ia ! <-- ! type of the boundary faces                     !
! ifrlag(nfabor)   ! ia ! --> ! type of the Lagrangian boundary faces          !
! itepa            ! ia ! <-- ! particle information (integers)                !
! (nbpmax,nivep    !    !     !                                                !
! xxpart           !  r ! <-- ! x-coordinate of the particle                   !
! yypart           !  r ! <-- ! y-coordinate of the particle                   !
! zzpart           !  r ! <-- ! z-coordinate of the particle                   !
! tvpart           !  r ! <-- ! value of the volume fraction                   !
! uupart           !  r ! <-- ! x-component of particle velocity               !
! vvpart           !  r ! <-- ! y-component of particle velocity               !
! wwpart           !  r ! <-- ! z-component of particle velocity               !
! ddpart           !  r ! <-- ! particle diameter                              !
! ttpart           !  r ! <-- ! particle temperature                           !
! dt(ncelet)       ! ra ! <-- ! time step (per cell)                           !
! rtpa             ! ra ! <-- ! transported variables at the previous          !
! (ncelet,*)       !    !     ! time step                                      !
! propce           ! ra ! <-- ! physical properties at cell centers            !
! (ncelet,*)       !    !     !                                                !
! ettp             ! ra ! <-- ! array of the variables associated to           !
!  (nbpmax,nvp)    !    !     ! the particles at the current time step         !
! tepa             ! ra ! <-- ! particle information (real) (statis. weight..) !
! (nbpmax,nvep)    !    !     !                                                !
!__________________!____!_____!________________________________________________!

!     Type: i (integer), r (real), s (string), a (array), l (logical),
!           and composite types (ex: ra real array)
!     mode: <-- input, --> output, <-> modifies data, --- work array

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

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

use paramx
use numvar
use optcal
use cstnum
use cstphy
use entsor
use lagpar
use lagran
use ppppar
use ppthch
use cpincl
use mesh

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

implicit none

! Arguments


integer          idvar  , iepart , izone  , iclass

integer          nvar   , nscal
integer          nbpmax , nvp    , nvp1   , nvep  , nivep
integer          ntersl , nvlsta , nvisbr

integer          itypfb(nfabor) , itrifb(nfabor)
integer          itepa(nbpmax,nivep) , ifrlag(nfabor)

double precision xxpart , yypart , zzpart
double precision tvpart , uupart , vvpart , wwpart
double precision ddpart , ttpart

double precision dt(ncelet) , rtpa(ncelet,*)
double precision propce(ncelet,*)
double precision ettp(nbpmax,nvp) , tepa(nbpmax,nvep)

! Local variables


double precision pis6

! Local variables
!Duc variable
integer it, ii, ii2
double precision :: r, v2, x, y, v3
!real :: ro1, ro2, m_sg, U1_const, U2_const,alp2
double precision, dimension (1 : 15) :: ri2, vi2
double precision, dimension (1 : 15) :: ri3, vi3
!Outlet variables



!===============================================================================
ri2(1:15) = (/0.284d-3, 0.71d-3, 1.278d-3, 1.846d-3, 2.272d-3, 2.84d-3, 3.266d-3,3.834d-3,&
 4.26d-3, 4.828d-3, 5.396d-3, 5.822d-3, 6.39d-3, 6.816d-3, 7.1d-3/) 
vi2(1:15) = (/11.274d0, 11.301d0, 11.017d0, 11.233d0, 11.079d0, 11.060d0, 10.960d0, 10.725d0,&
 10.45d0, 10.065d0, 9.59d0, 9.072d0, 8.054d0, 6.557d0, 4.746d0/)
ri3(1:15) = (/0.0d0,5.68d-4, 9.94d-4, 1.562d-3, 1.988d-3, 2.556d-3, 2.982d-3, 3.55d-3, &
 4.118d-3, 4.544d-3, 5.112d-3, 5.538d-3, 6.106d-3, 6.532d-3, 6.816d-3/) 
vi3(1:15) = (/10.03d0, 10.061d0, 10.027d0, 10.101d0, 10.105d0, 9.995d0, 10.137d0, &
 9.837d0, 9.942d0, 9.786d0, 9.529d0, 9.274d0, 9.015d0, 8.628d0, 8.169d0/)
!===============================================================================
! 1. Memory management
!===============================================================================


!===============================================================================
! 2. Initialization
!===============================================================================


!===============================================================================
! 3. Profile for the volume fraction
!==============================================================================
!Propriété de la première classe de particule

if (iclass .eq. 1) then 

   if (idvar .eq. 0) then
     tvpart = 0.000005
   endif

!===============================================================================
! 4. Velocity profile
!===============================================================================
    if (idvar .eq. 1) then
      x = xxpart
      y = yypart
      r=sqrt(x*x+y*y)   
      it = 0     
      ii2 = 15
      do ii = 1, ii2
        if (r >= ri2(ii) .and. r < ri2(ii+1)) then
          it = ii
          v2 = vi2(it) + (r-ri2(it))*(vi2(it+1)-vi2(it))/(ri2(it+1)-ri2(it))
          !write(NFECRA,*) ' X2=', x, ' Y2=', y,' r2=',r,' W2=', v2       
        endif                       
      enddo
      uupart = 0.d0
      vvpart = 0.d0
      wwpart = v2     
    endif

!===============================================================================
! 5. Diameter profile
!===============================================================================

    if (idvar .eq. 2) then

      ddpart = 25.d-6

    endif


!===============================================================================
! 6. Temperature profile
!===============================================================================

    if (idvar .eq. 3) then

      ttpart = 20.d0

    endif

endif

!Propriété de la deuxième classe de particule

if (iclass .eq. 2) then


   if (idvar .eq. 0) then
     tvpart = 0.0005d0
   endif

!===============================================================================
! 4. Velocity profile
!===============================================================================
    if (idvar .eq. 1) then
      x = xxpart
      y = yypart
      r=sqrt(x*x+y*y)   
      it = 0     
      ii2 = 15
      do ii = 1, ii2
        if (r >= ri3(ii) .and. r < ri3(ii+1)) then
          it = ii
          v3 = vi3(it) + (r-ri3(it))*(vi3(it+1)-vi3(it))/(ri3(it+1)-ri3(it))
          !write(NFECRA,*) ' X3=', x, ' Y3=', y,' r3=',r,' W3=', v3       
        endif                       
      enddo
      uupart = 0.d0
      vvpart = 0.d0
      wwpart = v3      
    endif

!===============================================================================
! 5. Diameter profile
!===============================================================================

    if (idvar .eq. 2) then

      ddpart = 70.d-6

    endif


!===============================================================================
! 6. Temperature profile
!===============================================================================

    if (idvar .eq. 3) then

      ttpart = 20.d0

    endif

endif

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

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

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

return

end subroutine uslapr


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

