 SUBROUTINE gradlim(msh)
!=========================== DOCUMENTATION ============================= 
!DOCDEBh*    GEO/gradlim.f90  *
!
!    NOM
! gradlim(msh)
!    DESCRIPTION  
!*
! Limite le gradient des variables primitives par la methode de Barth
! 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 a MPI pour mettre a jours TOUS les limiteurs
!
! Limits the gradient of primitive variables using Barth s method
! Therefore assumes that the gradient has already been calculated (call gradprim)
! Warning: the gradient of 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 :: dlim ! tableau des limiteurs de pente dans chaque cellule
  integer                                          :: i,iv,ic1,ic2,idom1,idom2,idom,ii,j
  real(kind=kind(0.d0)),dimension(nvarmax)         :: odg
  real(kind=kind(0.d0))                            :: g,r,dlim2
  type(vect_3)                                     :: vcf
  type(state)                                      :: w1
  type(state)                                      :: w2

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

  allocate(dlim(3,nvar,msh%nb_cell)) ! valeur min, valeur max, limiteur

  do i=1, msh%nb_cell
    do j=1, nvar
       dlim(1,j,i)=1.d20  !valeur min sur la cellule i de la variable j
       dlim(2,j,i)=-1.d20 !valeur max sur la cellule i de la variable j
       dlim(3,j,i)=1.d0   !limiteur  de la variable j de la cellule i
    enddo
  enddo

! 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
!-------------------------------------------------------------------------
! Determination du min et du max entourant une cellule
!-------------------------------------------------------------------------
! boucle sur les faces du maillage
  DO i=1,msh%nb_face
    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
    DO iv=1,nvar
      if(w2%vprim%v(iv).lt.dlim(1,iv,ic1))dlim(1,iv,ic1)=w2%vprim%v(iv)
      if(w2%vprim%v(iv).gt.dlim(2,iv,ic1))dlim(2,iv,ic1)=w2%vprim%v(iv)         
      IF (ic2.gt.0)then
        if(w1%vprim%v(iv).lt.dlim(1,iv,ic2))dlim(1,iv,ic2)=w1%vprim%v(iv)
        if(w1%vprim%v(iv).gt.dlim(2,iv,ic2))dlim(2,iv,ic2)=w1%vprim%v(iv)
      endif         
    END DO
  END DO
  DO i=1,msh%nb_cell
    w1=msh%list_cell(i)%w
    DO iv=1,nvar
      if(w1%vprim%v(iv).lt.dlim(1,iv,i))dlim(1,iv,i)=w1%vprim%v(iv)
      if(w1%vprim%v(iv).gt.dlim(2,iv,i))dlim(2,iv,i)=w1%vprim%v(iv)         
    END DO
  END DO

!-------------------------------------------------------------------------
! calcul du limiteur de pente de Barth-Jespersen
!-------------------------------------------------------------------------
! boucle sur les faces du maillage
  DO i=1,msh%nb_face
      
      ! reperage 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
         
         ! determination 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
             r=0.d0
             IF (g.ne.0.d0) THEN   ! moins compliqué et plus stable
                 if(g.gt.0.d0)r=(dlim(2,iv,ic1)-w1%vprim%v(iv))/g
                 if(g.lt.0.d0)r=(dlim(1,iv,ic1)-w1%vprim%v(iv))/g
                 IF (r.lt.0.d0) THEN
                     dlim2=0.d0
                     !write(*,*)'r negatif',ic1,r
                 ELSE IF (r.gt.1.d0) THEN
                     dlim2=1.d0
                     !write(*,*)'r sup a 1',ic1,r
                 ELSE
                     dlim2=r
                 END IF
             ELSE
                 dlim2=1.d0
             END IF
         
             IF (dlim(3,iv,ic1).gt.dlim2) dlim(3,iv,ic1)=dlim2
         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)
                 r=0.d0            
                 IF (g.ne.0.d0) THEN   ! moins compliqué et plus stable
 !                   if (dabs(g/odg(iv)).gt.1.d-6) then   ! si trop limité
                 if(g.ge.0.d0)r=(dlim(2,iv,ic2)-w2%vprim%v(iv))/g
                 if(g.lt.0.d0)r=(dlim(1,iv,ic2)-w2%vprim%v(iv))/g
                     IF (r.lt.0.d0) THEN
                         dlim2=0.d0
                         !write(*,*)'r negatif',ic2,r
                     ELSE IF (r.gt.1.d0) THEN
                         dlim2=1.d0
                         !write(*,*)'r sup a 1',ic2,r
                     ELSE
                         dlim2=r
                     END IF
                 ELSE
                  dlim2=1.d0
                 END IF
               
                 IF (dlim(3,iv,ic2).gt.dlim2) dlim(3,iv,ic2)=dlim2
               
             END DO
         END IF
      END IF
  END DO


! Mise a 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
              DO ii=1,msh%ndim
                  msh%list_cell(i)%w%gprim(ii)%v(iv)=0.9d0*dlim(3,iv,i)*msh%list_cell(i)%w%gprim(ii)%v(iv) ! précaution en plus...
                  if(msh%list_cell(i)%w%chi.ne.0.d0)msh%list_cell(i)%w%gprim(ii)%v(iv)=0.d0
              END DO
          END DO
      END IF
  END DO
            
  deallocate(dlim)
  
  
!===========================   FIN DE LA ROUTINE    ====================
END SUBROUTINE GRADLIM
!DOCFIN