 SUBROUTINE calc_flux(calc,lev)
!=========================== DOCUMENTATION ============================= 
!DOCDEBh*    NUM/calc_flux  *
!
!    NOM
!   calc_flux(calc,lev)
!    DESCRIPTION  
!*
! Calcul des bilan des flux sur les faces 
! Le terme de bilan est stocke dans la variable "balance" de chaque cellule
!
! Calculation of flow balances on faces 
! The balance term is stored in the ‘balance’ variable in each cell.
!*
! 
!    ENTREES / INPUT
!   calc : objet calcul / Calcul object
!   lev  :  niveau maximal de raffinement a considerer / maximum level of refinement to be considered
!    SORTIES / OUTPUT
!   calc : objet calcul avec les bilans de flux mis a jour / Calcul object with updated flow balances
!   lev  : niveau maximal de raffinement a considerer / maximum level of refinement to be considered
!
!    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 num_typ
  IMPLICIT NONE
!.2-----  Declaration
  type(calcul),intent(inout)  :: calc !! Objet calcul / Calcul object
  integer, intent(inout)      :: lev  !! niveau maximal de raffinement a considerer / maximum level of refinement to be considered
  
  type(state)                 :: wL,wR
  integer                     :: i,ic1,ic2,iv,ii
  real(kind=kind(0.d0))       :: surf,xa,ya,za
  type(vect_nvar)             :: flux, fluxnc 
  type(vect_3)                :: vnmoins ! oppose du vecteur normal (pour les flux non conservatifs)
  type(vect_3)                :: vcf     ! vecteur centre cellule - centre face
!=========================== DEBUT DU CODE EXECUTABLE ==================
  

!--- boucle sur les faces pour calculer les flux convectifs
  
  DO i=1,calc%msh%nb_face
      IF (calc%msh%list_face(i)%flevel.gt.lev) GOTO 200

      ic1=calc%msh%list_face(i)%ic1
      ic2=calc%msh%list_face(i)%ic2

      IF (ic2.eq.-999) GOTO 200
      surf=calc%msh%list_face(i)%surf
      IF (ic1.le.0) STOP
  
      ! recuperation de l etat Left
      wL=calc%msh%list_cell(ic1)%w
  
      ! si ordre 2, correction de l etat Left par la methode MUSCL
      IF (calc%iordre_e.ge.2) THEN
          ! determination du vecteur centre de la cellule - centre de la face
          xa=calc%msh%list_face(i)%center%x
          ya=calc%msh%list_face(i)%center%y
          za=calc%msh%list_face(i)%center%z
  
          vcf%v(1)=xa-calc%msh%list_cell(ic1)%center%x
          vcf%v(2)=ya-calc%msh%list_cell(ic1)%center%y
          vcf%v(3)=za-calc%msh%list_cell(ic1)%center%z
      
          DO iv=1,nvarmax
              DO ii=1,calc%msh%ndim
                  wL%vprim%v(iv)=wL%vprim%v(iv)+wL%gprim(ii)%v(iv)*vcf%v(ii) 
              END DO
          END DO
      
      END IF
   
      ! si il existe, recuperation de l etat Right
      IF (ic2.ge.1) THEN
          wR=calc%msh%list_cell(ic2)%w
           
         ! si ordre 2, correction de l etat Right par la methode MUSCL
          IF (calc%iordre_e.ge.2) THEN
              ! determination du vecteur centre de la cellule - centre de la face
              xa=calc%msh%list_face(i)%center%x
              ya=calc%msh%list_face(i)%center%y
              za=calc%msh%list_face(i)%center%z
        
              vcf%v(1)=xa-calc%msh%list_cell(ic2)%center%x
              vcf%v(2)=ya-calc%msh%list_cell(ic2)%center%y
              vcf%v(3)=za-calc%msh%list_cell(ic2)%center%z
              
              DO iv=1,nvarmax
                  DO ii=1,calc%msh%ndim
                    wR%vprim%v(iv)=wR%vprim%v(iv)+wR%gprim(ii)%v(iv)*vcf%v(ii)  
                  END DO
              END DO             
          END IF
      ELSE
          ! en cas de condition aux limites on extrapole l etat
          ! Right a partir du type de la condition aux limites
          ! de l etat connu et du vecteur normal
          xa=calc%msh%list_face(i)%center%x
          ya=calc%msh%list_face(i)%center%y
          za=calc%msh%list_face(i)%center%z
          CALL wext(wL,wR,calc%msh%list_face(i)%vnorm,ic2,xa,ya,za,&
                    calc%msh%list_cell(ic1)%center%x, &
                    calc%msh%list_cell(ic1)%center%y, &
                    calc%msh%list_cell(ic1)%center%z )
      END IF
  
      ! calcul du flux numerique convectif conservatif
      icell_phy=ic1

      CALL prim2bal(wR)
      CALL prim2bal(wL)
      CALL fluxnum(wL,wR,calc%msh%list_face(i)%vnorm,flux,calc%ityp_flux)

      ! ajout dans le bilan de la cellule de gauche si elle est dans le
      ! bon domaine
      IF (calc%msh%list_cell(ic1)%idom.eq.calc%msh%numdom.and.calc%msh%list_cell(ic1)%isolid.ne.1) &
             calc%msh%list_cell(ic1)%balance%v=calc%msh%list_cell(ic1)%balance%v+flux%v * surf
                   
      ! ajout dans le bilan de la cellule de droite si elle existe et est dans le
      ! bon domaine (signe oppose car f(-vn)=-f(vn) )
      IF (ic2.gt.0) THEN
          IF (calc%msh%list_cell(ic2)%idom.eq.calc%msh%numdom.and.calc%msh%list_cell(ic2)%isolid.ne.1) &
                 calc%msh%list_cell(ic2)%balance%v=calc%msh%list_cell(ic2)%balance%v-flux%v *surf
      END IF
      ! le cas echant, calcul du flux numerique non conservatif

          CALL fluxnumnoncons(wL,wR,calc%msh%list_face(i)%vnorm,fluxnc,calc%ityp_flux)
          
          ! ajout dans le bilan de la cellule de gauche si elle est dans le
          ! bon domaine
          IF (calc%msh%list_cell(ic1)%idom.eq.calc%msh%numdom) &
              calc%msh%list_cell(ic1)%balance%v=calc%msh%list_cell(ic1)%balance%v+fluxnc%v * surf
  
          vnmoins%v=-calc%msh%list_face(i)%vnorm%v
          vi_phy=-vi_phy
          call fluxnumnoncons(wR,wL,vnmoins,fluxnc,calc%ityp_flux)
          ! ajout dans le bilan de la cellule de droite si elle existe et est dans le
          ! bon domaine (ici en general, f(-vn) est different de -f(vn))
          IF (ic2.gt.0) THEN
             IF (calc%msh%list_cell(ic2)%idom.eq.calc%msh%numdom) &
                  calc%msh%list_cell(ic2)%balance%v=calc%msh%list_cell(ic2)%balance%v+fluxnc%v *surf
          END IF

200 continue
  END DO
  
!===========================   FIN DE LA ROUTINE    ====================
END SUBROUTINE calc_flux
!DOCFIN