 SUBROUTINE Wext(win,wout,vn,ityp,x,y,z,xm,ym,zm)
!=========================== DOCUMENTATION ============================= 
!DOCDEBh*    PHY/Wext  *
!
!    NOM
! wext(win,wout,vn,ityp)
!    DESCRIPTION  
!*
! Determination des conditions aux limites
! Determining boundary conditions
!*
! 
!    ENTREES / INPUT
! win  : objet de type "state" contenant les variables de la cellule courante 
!      / state object containing the variables of the current cell
! vn   : objet de type "vect_3", normale a la face / normal to the face
! ityp : Type de condition aux limites / Type of boundary condition
!        0 -> On recopie / We copy
!       -1 -> Condition miroir / Mirror condition
! x,y,z : centre de la facette ou la condition est calculee / center of the facet where the condition is calculated
! xm,ym,zm : centre de la cellule de win. / center of the cell of win
!
!    SORTIES / OUTPUT
! wout : objet de type "state" contenant la condition aux limites
!
!    SOURCE
! Ce logiciel est regi par la licence [CeCILL-B]
! This software is governed by the [CeCILL-B] license
!=========================== DEBUT DES DECLARATIONS ==================== 
!.1-----  Implicit, Use
  use phy_typ
  use phy ,only: prim2bal, fct_cond,phy_state,bal2prim
  IMPLICIT NONE
!.2-----  Declaration
  integer,intent(in)                        :: ityp     !! Type de condition aux limites / Type of boundary condition
  type(state),intent(in)                    :: win      !! objet de type "state" de la cellule courante / state object of the current cell
  type(vect_3),intent(in)                   :: vn       !! normale a la face / normal to the face
  type(state),intent(out)                   :: wout     !! objet de type "state" sur la cellule fictive / state object on the fictitious cell
  real(kind=kind(0.d0)),optional,intent(in) :: x,y,z    !! position du point sur la frontiere / position of point on boundary
  real(kind=kind(0.d0)),optional,intent(in) :: xm,ym,zm !! centre de la cellule courante / centre of current cell 
  
  real(kind=kind(0.d0))                     :: qn,q
  integer                                   :: fctnum,i
  !=========================== DEBUT DU CODE EXECUTABLE ==================

wout%vprim%v=0.d0
wout%vbal%v=0.d0
wout%gprim(1)%v=0.d0
wout%gprim(2)%v=0.d0
wout%gprim(3)%v=0.d0

  SELECT CASE (ityp)
      CASE(0) !  recopie simple
          wout=win

      CASE (-1) ! condition de non-penetration
          wout%vprim=win%vprim
          qn=win%vprim%v(2)*vn%v(1)+win%vprim%v(3)*vn%v(2)+win%vprim%v(4)*vn%v(3)
          wout%vprim%v(2)=win%vprim%v(2)-2.d0*qn*vn%v(1)
          wout%vprim%v(3)=win%vprim%v(3)-2.d0*qn*vn%v(2)
          wout%vprim%v(4)=win%vprim%v(4)-2.d0*qn*vn%v(3)            
          wout%vprim%v(5)=win%vprim%v(5)-2.d0*win%vprim%v(1)*gpes_phy*(y-ym)
          if(imodel_phy.eq.1)then
              wout%vprim%v(4)=0.d0*xm*zm*ym       
              wout%vprim%v(5)=win%vprim%v(5)
          endif
          CALL prim2bal(wout)

      CASE (-2) ! condition d adherence
          wout%vprim=win%vprim
          wout%vprim%v(2)=0.d0
          wout%vprim%v(3)=0.d0
          wout%vprim%v(4)= 0.d0           
          CALL prim2bal(wout)  

        
! conditions imposees: on cherche l etat impose dans le vecteur phy_state

      CASE (-101)
          wout=phy_state(1)  
          !on laisse libre les ddls dont phy_state est > 1.d20
          DO i=1,nvar
              IF (phy_state(1)%vprim%v(i).ge.1.d20) wout%vprim%v(i)= win%vprim%v(i)
          END DO
          if(imodel_phy.eq.1)wout%vprim%v(5)=win%vprim%v(5)
          CALL prim2bal(wout)

      CASE (-102)
          wout=phy_state(2) 
          !on laisse libre les ddls dont phy_state est > 1.d20
          DO i=1,nvar
              IF (phy_state(2)%vprim%v(i).ge.1.d20) wout%vprim%v(i)= win%vprim%v(i)
          END DO
          if(imodel_phy.eq.1)wout%vprim%v(5)=win%vprim%v(5)
          CALL prim2bal(wout)

      CASE (-103)
          wout=phy_state(3) 
          !on laisse libre les ddls dont phy_state est > 1.d20
          DO i=1,nvar
              IF (phy_state(3)%vprim%v(i).ge.1.d20) wout%vprim%v(i)= win%vprim%v(i)
          END DO
          if(imodel_phy.eq.1)wout%vprim%v(5)=win%vprim%v(5)
          CALL prim2bal(wout)

      CASE (-104)
          wout=phy_state(4) 
          !on laisse libre les ddls dont phy_state est > 1.d20
          DO i=1,nvar
              IF (phy_state(4)%vprim%v(i).ge.1.d20) wout%vprim%v(i)= win%vprim%v(i)
          END DO
          if(imodel_phy.eq.1)wout%vprim%v(5)=win%vprim%v(5)
          CALL prim2bal(wout)

      CASE (-105)
          wout=phy_state(5) 
          !on laisse libre les ddls dont phy_state est > 1.d20
          DO i=1,nvar
              IF (phy_state(5)%vprim%v(i).ge.1.d20) wout%vprim%v(i)= win%vprim%v(i)
          END DO
          CALL prim2bal(wout)
      
      CASE (-106)
          wout=phy_state(6) 
          !on laisse libre les ddls dont phy_state est > 1.d20
          DO i=1,nvar
              IF (phy_state(6)%vprim%v(i).ge.1.d20) wout%vprim%v(i)= win%vprim%v(i)
          END DO
          CALL prim2bal(wout)
      
      CASE (-107)
          wout=phy_state(7) 
      
      CASE (-108)
          wout=phy_state(8) 
      
      CASE (-109)
          wout=phy_state(9) 
      
      CASE (-110)
          wout=phy_state(10) 
          
      CASE (-201)
          wout=phy_state(1)  
          !on laisse libre les ddls dont phy_state est > 1.d20
          DO i=1,nvar
              IF (phy_state(1)%vbal%v(i).ge.1.d20) wout%vbal%v(i)= win%vbal%v(i)
          END DO
          if(imodel_phy.eq.1)wout%vprim%v(5)=win%vprim%v(5)
          q=sqrt(phy_state(1)%vbal%v(2)**2+phy_state(1)%vbal%v(3)**2)
         if(imodel_phy.eq.1.and.win%vbal%v(1).le.hzero_phy.and.q.gt.1.d-6) then
           wout%vbal%v(1)=(q**2/gpes_phy)**(1.d0/3.d0)
         endif
          CALL bal2prim(wout)

      CASE (-202)
          wout=phy_state(2) 
          !on laisse libre les ddls dont phy_state est > 1.d20
          DO i=1,nvar
              IF (phy_state(2)%vbal%v(i).ge.1.d20) wout%vbal%v(i)= win%vbal%v(i)
          END DO
          if(imodel_phy.eq.1)wout%vprim%v(5)=win%vprim%v(5)
          CALL bal2prim(wout)

      CASE (-203)
          wout=phy_state(3) 
          !on laisse libre les ddls dont phy_state est > 1.d20
          DO i=1,nvar
              IF (phy_state(3)%vbal%v(i).ge.1.d20) wout%vbal%v(i)= win%vbal%v(i)
          END DO
          if(imodel_phy.eq.1)wout%vprim%v(5)=win%vprim%v(5)
          CALL bal2prim(wout)

      CASE (-204)
          wout=phy_state(4) 
          !on laisse libre les ddls dont phy_state est > 1.d20
          DO i=1,nvar
              IF (phy_state(4)%vbal%v(i).ge.1.d20) wout%vbal%v(i)= win%vbal%v(i)
          END DO
          if(imodel_phy.eq.1)wout%vprim%v(5)=win%vprim%v(5)
          CALL bal2prim(wout)

      CASE (-205)
          wout=phy_state(5) 
          !on laisse libre les ddls dont phy_state est > 1.d20
          DO i=1,nvar
              IF (phy_state(5)%vbal%v(i).ge.1.d20) wout%vbal%v(i)= win%vbal%v(i)
          END DO
          CALL bal2prim(wout)
      
      CASE (-206)
          wout=phy_state(6) 
          !on laisse libre les ddls dont phy_state est > 1.d20
          DO i=1,nvar
              IF (phy_state(6)%vbal%v(i).ge.1.d20) wout%vbal%v(i)= win%vbal%v(i)
          END DO
          CALL bal2prim(wout)
      


      CASE (-999) ! frontiere artificielle pour le passage du 3D au 2D, on fait comme un miroir            
          wout%vprim=win%vprim
          qn=win%vprim%v(2)*vn%v(1)+win%vprim%v(3)*vn%v(2)+win%vprim%v(4)*vn%v(3)
          wout%vprim%v(2)=win%vprim%v(2)-2.d0*qn*vn%v(1)
          wout%vprim%v(3)=win%vprim%v(3)-2.d0*qn*vn%v(2)
          wout%vprim%v(4)=win%vprim%v(4)-2.d0*qn*vn%v(3)            
          CALL prim2bal(wout)

      CASE (-998) ! frontiere artificielle pour le passage du 3D au AXIsymetrique, on fait tourner la quantite de mouvement             
          wout=win
          wout%vprim%v(2)=win%vprim%v(2)*(vn%v(3)**2-vn%v(1)**2)
          wout%vprim%v(4)=-2.d0*win%vprim%v(2)*vn%v(3)*vn%v(1)           
          CALL prim2bal(wout)

      CASE (-20000:-10000) ! Cas ou la condition aux limites est definie par une fonction utilisateur
          wout=win
          fctnum=-10000-ityp
          CALL fct_cond(win%vprim,wout%vprim,fctnum,x,y,z,vn)
          CALL prim2bal(wout)
           
      CASE DEFAULT
          write(*,*)ityp 
          call print_err('WEXT','Kind of boundary condition not defined')
        
  END SELECT


!===========================   FIN DE LA ROUTINE    ====================
END SUBROUTINE Wext
!DOCFIN