 SUBROUTINE gradlimxy(msh)
!=========================== DOCUMENTATION ============================= 
!DOCDEBh*    GEO/gradlimxy.f90  *
!
!    NOM
! gradlimxy(msh)
!    DESCRIPTION  
!*
! Limite le gradient des variables primitives par la méthode de Barth
! La limitation est faite dans chaque direction x y indépendamment
! Suppose donc que le maillage est régulier en x y et 2D 
! Suppose donc que le gradient  a deja ete calcule (appel a gradprim)
! Attention: le gradient des cellules frontiere n est pas limite
! Il faut donc un appel à MPI pour mettre à jours TOUS les limiteurs
!
! Limit the gradient of primitive variables using Barth s method
! Limitation is done in each x y direction independently
! Assumes that the mesh is regular in x y and 2D. 
! Assumes that the gradient has already been calculated (call gradprim).
! Warning: the gradient of the border cells is not limited
! A call to MPI is therefore required to update ALL the limiters.
!*
! 
!    ENTREES / INPUT
!    msh : maillage / mesh
!
!    SORTIES / OUTPUT
!    msh : maillage / mesh
!
!    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 geo_typ
  IMPLICIT NONE
!.2-----  Declaration
  type(mesh), intent(inout) :: msh !! Maillage / Mesh
  real(kind=kind(0.d0)),dimension(:,:),allocatable :: dlimx , dlimy ! tableau des limiteurs de pente dans chaque cellule
  real(kind=kind(0.d0)),dimension(nvarmax)         :: odg
  integer                                          :: i,iv,ic1,ic2,idom1,idom2,idom
  real(kind=kind(0.d0))                            :: g,r,dlim2
  type(vect_3)                                     :: vcf
  type(state)                                      :: w1,w2

!=========================== DEBUT DU CODE EXECUTABLE ==================

  allocate(dlimx(nvar,msh%nb_cell))
  allocate(dlimy(nvar,msh%nb_cell))

  dlimx=1
  dlimy=1

! vecteur des ordres de grandeur (permet d avoir un limiteur adimensionne)
  DO iv=1,nvar
      odg(iv)=1.d0
  END DO

  odg(5)=1.d5
!  odg(1)=1.d0
  
  ! boucle sur les faces du maillage
  DO i=1,msh%nb_face
      
      ! repérage des voisins
      ic1=msh%list_face(i)%ic1
      ic2=msh%list_face(i)%ic2
      
      w1=msh%list_cell(ic1)%w
      IF (ic2.gt.0) THEN
          w2=msh%list_cell(ic2)%w
      ELSE
          CALL wext(w1,w2,msh%list_face(i)%vnorm,ic2, &
               msh%list_face(i)%center%x , &
               msh%list_face(i)%center%y , &
               msh%list_face(i)%center%z , &
               msh%list_cell(ic1)%center%x, &
               msh%list_cell(ic1)%center%y, &
               msh%list_cell(ic1)%center%z )
      END IF
     
      ! numero de domaine à gauche
      idom1=msh%list_cell(ic1)%idom
      IF (idom1.ne.msh%numdom) THEN
          write(*,*) 'erreur maillage domaine'
          STOP
      END IF
     
      IF (ic2.ne.-999) THEN ! ne pas boucler sur les faces fictives
        
          ! numero de domaine à droite
          idom2=0
          IF (ic2.gt.0) idom2=msh%list_cell(ic2)%idom
          IF (idom2.ne.msh%numdom) idom2=0
          
          
          ! limitation cellule de gauche
          
          ! détermination du vecteur centre de la cellule - centre de la face
          vcf%v(1)=msh%list_face(i)%center%x-msh%list_cell(ic1)%center%x
          vcf%v(2)=msh%list_face(i)%center%y-msh%list_cell(ic1)%center%y
          vcf%v(3)=msh%list_face(i)%center%z-msh%list_cell(ic1)%center%z
       
          DO iv=1,nvar
              
              dlim2=1.d0
              g=w1%gprim(1)%v(iv)*vcf%v(1)+&
              & w1%gprim(2)%v(iv)*vcf%v(2)+&
              & w1%gprim(3)%v(iv)*vcf%v(3)
              
              !IF (dabs(g/odg(iv)).gt.1.d-6) THEN
              IF (g.ne.0.d0) THEN   ! moins compliqué et plus stable
                  r=(w2%vprim%v(iv)-w1%vprim%v(iv))/g
                  IF (r.lt.0.d0) THEN
                      dlim2=0.d0
                  ELSE IF (r.gt.1.d0) THEN
                      dlim2=1.d0
                  ELSE
                      dlim2=r
                  END IF
              ELSE
              dlim2=1.d0
              END IF
              
              IF (dabs(msh%list_face(i)%vnorm%v(2)).lt.1.d-10) THEN         ! direction x
                  IF (dlimx(iv,ic1).gt.dlim2) dlimx(iv,ic1)=dlim2
              END IF
              
              IF (dabs(msh%list_face(i)%vnorm%v(1)).lt.1.d-10) THEN         ! direction y
                  IF (dlimy(iv,ic1).gt.dlim2) dlimy(iv,ic1)=dlim2
              END IF
              
          END DO
        
          ! limitation cellule de droite
        
          IF (idom2.ne.0) THEN
           
              ! détermination du vecteur centre de la cellule - centre de la face
              vcf%v(1)=msh%list_face(i)%center%x-msh%list_cell(ic2)%center%x
              vcf%v(2)=msh%list_face(i)%center%y-msh%list_cell(ic2)%center%y
              vcf%v(3)=msh%list_face(i)%center%z-msh%list_cell(ic2)%center%z
               
              DO iv=1,nvar
                  
                  dlim2=1.d0
                  
                  g=w2%gprim(1)%v(iv)*vcf%v(1)+&
                  & w2%gprim(2)%v(iv)*vcf%v(2)+&
                  & w2%gprim(3)%v(iv)*vcf%v(3)
                  
                  IF (g.ne.0.d0) THEN   ! moins compliqué et plus stable
!                  if (dabs(g/odg(iv)).gt.1.d-6) then   ! si trop limité
                      r=(w1%vprim%v(iv)-w2%vprim%v(iv))/g
                      IF (r.lt.0.d0) THEN
                          dlim2=0.d0
                      ELSE IF (r.gt.1.d0) THEN
                          dlim2=1.d0
                      ELSE
                          dlim2=r
                      END IF
                  ELSE
                     dlim2=1.d0
                  END IF
                  
                  IF (dabs(msh%list_face(i)%vnorm%v(2)).lt.1.d-10) THEN         ! direction x
                     IF (dlimx(iv,ic2).gt.dlim2) dlimx(iv,ic2)=dlim2
                  END IF
                  
                  IF (dabs(msh%list_face(i)%vnorm%v(1)).lt.1.d-10) THEN         ! direction y
                     IF (dlimy(iv,ic2).gt.dlim2) dlimy(iv,ic2)=dlim2
                  END IF
                  
              END DO

          END IF
        
      END IF
     
  END DO


  ! mise à jour du gradient dans les cellules

  DO i=1,msh%nb_cell
      idom=msh%list_cell(i)%idom
      IF (idom.eq.msh%numdom) THEN
          DO iv=1,nvar
              msh%list_cell(i)%w%gprim(1)%v(iv)=0.95*dlimx(iv,i)*msh%list_cell(i)%w%gprim(1)%v(iv) 
              msh%list_cell(i)%w%gprim(2)%v(iv)=0.95*dlimy(iv,i)*msh%list_cell(i)%w%gprim(2)%v(iv) 
          END DO
      END IF  
  END DO

  deallocate(dlimx)
  deallocate(dlimy)
!===========================   FIN DE LA ROUTINE    ====================
END SUBROUTINE gradlimxy
!DOCFIN