TABLE OF CONTENTS

GEO/gradprim [ Modules ]

[ Top ] [ Modules ]

NOM

 gradprim(msh)

DESCRIPTION

 Calcul du gradient des variables primitives sur les cellules du maillage msh
 dans le domaine courant.
 attention : les gradients des cellules frontieres ne sont pas calcules

 Calculation of the gradient of the primitive variables on the cells of the msh mesh
 in the current domain.
 Warning: the gradients of border cells are not calculated
 
    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
  use geo, only: display
  IMPLICIT NONE
!.2-----  Declaration
  type(mesh), intent(inout)    :: msh !! Maillage / Mesh
  
  type(state)                  :: w1,w2
  type(vect_nvar),dimension(3) :: vpn
  integer                      :: i,ii,ic1,ic2,idom1,idom2  
!=========================== DEBUT DU CODE EXECUTABLE ==================

  DO i=1,msh%nb_cell
      DO ii=1,3
          msh%list_cell(i)%w%gprim(ii)%v(1:nvar)=0.d0
      END DO
      
!SV: on stocke eta dans la variable z. Retour à la normale a la fin.
      IF(imodel_phy.eq.1) &
      msh%list_cell(i)%w%vprim%v(5)=msh%list_cell(i)%w%vprim%v(1)+msh%list_cell(i)%w%vprim%v(5)
  END DO

  DO i=1,msh%nb_face
      ic1=msh%list_face(i)%ic1
      ic2=msh%list_face(i)%ic2
      
      IF (ic1.le.0) THEN
          write(*,*) 'erreur maillage orientation'
          STOP
      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
          
          ! états gauche et droit
          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
          
          ! formule de green pour calculer le gradient au centre
          ! int(Omega,grad f)= int(dOmega, f*normale) = int(dOmega,1/2( f*norm(L) +f*norm(R)))
          ! ce calcul devra être modifié pour des maillages très irréguliers 
          ! ou si on passe a l ordre 2 !!!
          DO ii=1,msh%ndim
              vpn(ii)%v=0.5d0*msh%list_face(i)%vnorm%v(ii)*msh%list_face(i)%surf*&
                  &(w1%vprim%v+w2%vprim%v)
             
              msh%list_cell(ic1)%w%gprim(ii)%v(1:nvar)=msh%list_cell(ic1)%w%gprim(ii)%v(1:nvar)+&
                  &1.d0/msh%list_cell(ic1)%vol*vpn(ii)%v(1:nvar)
             if(imodel_phy.eq.1.and.ngrad_bathy.eq.1)msh%list_cell(i)%w%gprim(ii)%v(5) = msh%list_cell(ic1)%w%gprim(ii)%v(5)+&
                  &1.d0/msh%list_cell(ic1)%vol*vpn(ii)%v(5)
              IF (idom2.ne.0) THEN
                  msh%list_cell(ic2)%w%gprim(ii)%v(1:nvar)=msh%list_cell(ic2)%w%gprim(ii)%v(1:nvar)-&
                     &1.d0/msh%list_cell(ic2)%vol*vpn(ii)%v(1:nvar)
                  if(imodel_phy.eq.1.and.ngrad_bathy.eq.1)msh%list_cell(i)%w%gprim(ii)%v(5) = msh%list_cell(ic2)%w%gprim(ii)%v(5)+&
                       &1.d0/msh%list_cell(ic1)%vol*vpn(ii)%v(5)
              END IF

          END DO
      
      END IF
  
  END DO

 !===========================!=========================== !=========================== !=========================== 
! on calcule le gradient de phi en fonction de grap p et grap rho
!
  DO i=1,msh%nb_cell
!         do ii=1,3
!            if (imodel_phy.eq.3  .OR. imodel_phy.eq.8) then
!            msh%list_cell(i)%w%gprim(ii)%v(10) =-1.D0/((rho0W_phy-rho0A_phy)*cW_phy**2)&
!           *(msh%list_cell(i)%w%gprim(ii)%v(5)-cW_phy**2*msh%list_cell(i)%w%gprim(ii)%v(1))
!            else if (imodel_phy.eq.4) then
!            msh%list_cell(i)%w%gprim(ii)%v(10) =-1.D0/((rho0W_phy-rho0A_phy)*cW_phy**2)&
!           *(msh%list_cell(i)%w%gprim(ii)%v(6)-cW_phy**2*msh%list_cell(i)%w%gprim(ii)%v(1))
!            END IF
!          END DO
!SV : retour aux variables initiales (z=(z+h)-h et idem pour grad_z)
!     on tue les grad_z de Zbottom pour eviter perturber le calcul des limiteurs
!     (mais potentiellement non necessaire!! => certainement a supprimer apres tests!)
      IF (imodel_phy.eq.1)then
         msh%list_cell(i)%w%gprim(:)%v(4)=0.d0
! dans le cas ou le gradient de la bathymetrie est numerique
         if(ngrad_bathy.eq.1)msh%list_cell(i)%w%gprim(:)%v(5)=msh%list_cell(i)%w%gprim(:)%v(5)-& 
                                                              msh%list_cell(i)%w%gprim(:)%v(1)
         msh%list_cell(i)%w%vprim%v(5)=msh%list_cell(i)%w%vprim%v(5)-msh%list_cell(i)%w%vprim%v(1)
      ENDIF
  END DO
!
!===========================   FIN DE LA ROUTINE    ====================
    END SUBROUTINE gradprim