C YGMV      SOURCE    GOUNAND   25/11/12    21:15:50     12399          
      SUBROUTINE YGMV
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
C***********************************************************************
C
C  CE SP DISCRETISE LE TERME DE SOURCE DE QDM DANS LES EQUATIONS DE
C  NAVIER STOKES
C    EN 2D SUR LES ELEMENTS QUA4 ET TRI3      PLAN OU AXI
C    EN 3D SUR LES ELEMENTS CUB8 ET PRI6
C  L OPERATEUR EST "SOUS-INTEGRE"
C
C  SYNTAXE :
C  ---------
C
C         GMV TABGMV  INCO UN :
C
C  COEFFICIENT :
C  -------------
C
C  TABGMV   Table contenant les entrees suivantes
C
C
C  INCONNUES :
C  -----------
C
C   UN                  CHAMPS DE VITESSE
C
C
C***********************************************************************


-INC PPARAM
-INC CCOPTIO
-INC SMLREEL
      POINTEUR MQR.MLREEL,MPR.MLREEL
-INC SMEVOLL

-INC SMCOORD
-INC SMLENTI
      POINTEUR IZIPAD.MLENTI
-INC SMELEME
      POINTEUR MELEM1.MELEME,MELEMC.MELEME,MDEBI.MELEME
      POINTEUR MENTR.MELEME,MSORT.MELEME
-INC SMCHPOI
      POINTEUR IZG1.MCHPOI,IZGG1.MPOVAL
      POINTEUR IZTU1.MPOVAL,IZPP.MPOVAL
      POINTEUR IZVOL.MPOVAL


-INC SMTABLE
      POINTEUR KIZG.TABLE,MTABX.MTABLE,OPTI.MTABLE
      POINTEUR INCO.TABLE,KOPT.MTABLE
      POINTEUR MTABZ.MTABLE,MTABD.MTABLE,MTABP.MTABLE
      POINTEUR MTABA.MTABLE
-INC SMLMOTS
      POINTEUR LINCO.MLMOTS
      CHARACTER*8 NOMZ,NOMI,NOMA,TYPE,TYPC
      REAL*8 XVEC1(3),XVEC2(3)
      PARAMETER (NTB=1)
      CHARACTER*8 LTAB(NTB)
      DIMENSION KTAB(NTB)
      SAVE IPAS
      DATA LTAB/'KIZX    '/
      DATA IPAS/0/
C*****************************************************************************

CGMV
      CALL LITABS(LTAB,KTAB,NTB,1,IRET)
      IF(IRET.EQ.0)THEN
      WRITE(6,*)' Operateur GMV '
      WRITE(6,*)' On attend un ensemble de table soustypes'
      RETURN
      ENDIF
      MTABX=KTAB(1)
      SEGACT MTABX

C*****************************************************************************
C OPTIONS
C CES PARAMETRES SONT INITIALISES POUR ETRE EN DECENTRE ET NE PAS AVOIR DE
C POROSITE :

      IOP4=0
      IAXI=0
      IF(IFOMOD.EQ.0)IAXI=2
      CALL LEKTAB(MTABX,'KOPT',OPTI)
      IF(OPTI.NE.0)THEN
      SEGACT OPTI
      TYPE=' '
      CALL ACMO(OPTI,'POROSITE',TYPE,MCHPOI)
      IF(TYPE.EQ.'CHPOINT ')THEN
      IOP7=1
      CALL LICHT(MCHPOI,IZPORO,TYPC,IGEOM)
      ENDIF
      ENDIF
C*****************************************************************************

      CALL ACMM(MTABX,'NOMZONE',NOMZ)
      CALL LEKTAB(MTABX,'DOMZ',MTABZ)
      IF(MTABZ.EQ.0)THEN
      WRITE(6,*)' Operateur GMV '
      WRITE(6,*)' On ne trouve pas l''indice DOMZ ? '
      GO TO 90
      ENDIF
      SEGACT MTABZ

      CALL LEKTAB(MTABZ,'MAILLAGE',MELEME)
      IF(MELEME.EQ.0)GO TO 90
      SEGACT MELEME


      CALL LEKTAB(MTABZ,'XXVOLUM',MCHPOI)
      IF(MCHPOI.EQ.0)GO TO 90
      CALL LICHT(MCHPOI,IZVOL,TYPC,IGEOM)

C***

      TYPE='LISTMOTS'
      CALL ACMO(MTABX,'LISTINCO',TYPE,LINCO)
      SEGACT LINCO

      CALL ACME(MTABX,'IARG',IARG)

      IF(IARG.NE.1)THEN
      WRITE(6,*)' Operateur GMV '
      WRITE(6,*)' Nombre d''arguments ( ',IARG,' ) incorrect '
      RETURN
      ENDIF

      TYPE=' '
      CALL ACMO(MTABX,'ARG1',TYPE,MTABA)
      IF(TYPE.NE.'TABLE')THEN
      WRITE(6,*)' Operateur GMV '
      WRITE(6,*)
     &' LE TYPE DE L''ARGUMENT (',TYPE,') N EST PAS CONVENABLE'
      WRITE(6,*)' On attend une table '
      RETURN
      ELSE
      SEGACT MTABA

      TYPE=' '
      CALL ACMO(MTABA,'DIR',TYPE,IPQ)

      IF(TYPE.NE.'POINT')THEN
      WRITE(6,*)' Operateur GMV '
      WRITE(6,*)' Entree DIR erronee '
      RETURN
      ELSE
      XVEC1(1)=XCOOR((IPQ-1)*(IDIM+1)    +1)
      XVEC1(2)=XCOOR((IPQ-1)*(IDIM+1)    +2)
      XNN=XVEC1(1)*XVEC1(1)+XVEC1(2)*XVEC1(2)
      IF(IDIM.EQ.3)THEN
      XVEC1(3)=XCOOR((IPQ-1)*(IDIM+1)    +3)
      XNN=XNN+XVEC1(3)*XVEC1(3)
      ENDIF
      XVEC2(1)=XVEC1(1)/XNN
      XVEC2(2)=XVEC1(2)/XNN
      IF(IDIM.EQ.3)XVEC2(3)=XVEC1(3)/XNN

      ENDIF

      TYPE=' '
      CALL ACMO(MTABA,'IMPR',TYPE,IRET)
      IF(TYPE.NE.'ENTIER')THEN
      IMPR=0
      ELSE
      CALL ACME(MTABA,'IMPR',IMPR)
      ENDIF

      TYPE=' '
      CALL ACMO(MTABA,'KIMP',TYPE,IRET)
      IF(TYPE.NE.'FLOTTANT')THEN
      IKIMP=0
      ELSE
      IKIMP=1
      CALL ACMF(MTABA,'KIMP',AKIMP)
      ENDIF

      IF(IKIMP.EQ.0)THEN

        TYPE=' '
        CALL ACMO(MTABA,'K0',TYPE,IRET)
        IF(TYPE.NE.'FLOTTANT')THEN
        AK0=1.
        CALL ECMF(MTABA,'K0',AK0)
        ELSE
        CALL ACMF(MTABA,'K0',AK0)
        ENDIF

        TYPE=' '
        CALL ACMO(MTABA,'OMEGA',TYPE,IRET)
        IF(TYPE.NE.'FLOTTANT')THEN
        W1=0.1
        CALL ECMF(MTABA,'OMEGA',W1)
        ELSE
        CALL ACMF(MTABA,'OMEGA',W1)
        ENDIF

        TYPE=' '
        CALL ACMO(MTABA,'GMV',TYPE,MEVOLL)
        IF(TYPE.EQ.'EVOLUTIO')THEN
        SEGACT MEVOLL
        KEVOLL=IEVOLL(1)
        SEGACT KEVOLL
        MQR=IPROGX
        MPR=IPROGY
        SEGDES KEVOLL,MEVOLL
        ELSE
        WRITE(6,*)' Operateur GMV '
        WRITE(6,*)' Entree GMV erronee '
        WRITE(6,*)' On attend un type EVOLUTION'
        RETURN
        ENDIF

      ENDIF


      TYPE=' '
      CALL ACMO(MTABA,'LDEBIT',TYPE,MDEBI)
      IF(TYPE.NE.'MAILLAGE')THEN
      WRITE(6,*)' Operateur GMV '
      WRITE(6,*)' Entree LDEBIT erronee '
      WRITE(6,*)' On attend un type MAILLAGE'
      RETURN
      ENDIF

      TYPE=' '
      CALL ACMO(MTABA,'PENTREE',TYPE,MENTR)
      IF(TYPE.NE.'MAILLAGE')THEN
      WRITE(6,*)' Operateur GMV '
      WRITE(6,*)' Entree PENTREE erronee '
      WRITE(6,*)' On attend un type MAILLAGE'
      RETURN
      ENDIF

      TYPE=' '
      CALL ACMO(MTABA,'PSORTIE',TYPE,MSORT)
      IF(TYPE.NE.'MAILLAGE')THEN
      WRITE(6,*)' Operateur GMV '
      WRITE(6,*)' Entree PSORTIE erronee '
      WRITE(6,*)' On attend un type MAILLAGE'
      RETURN
      ENDIF

      ENDIF
C****************** Fin lecture Table Arguments *****************

      CALL LEKTAB(MTABX,'EQEX',MTAB1)
      IF(MTAB1.EQ.0)THEN
      WRITE(6,*)' Operateur GMV '
      WRITE(6,*)' On ne trouve pas l''indice EQEX ? '
      GO TO 90
      ENDIF
      SEGACT MTAB1
      CALL LEKTAB(MTAB1,'DOMAINE',MTABD)
      IF(MTABD.EQ.0)THEN
      WRITE(6,*)' Operateur GMV '
      WRITE(6,*)' On ne trouve pas l''indice DOMAINE ?'
      GO TO 90
      ENDIF
      CALL LEKTAB(MTABD,'SOMMET',MELEM1)
      IF(MELEM1.EQ.0)THEN
      WRITE(6,*)' Operateur GMV '
      WRITE(6,*)' On ne trouve pas l''indice SOMMET ?'
      GO TO 90
      ENDIF

      CALL LEKTAB(MTABD,'CENTRE',MELEMC)
      IF(MELEMC.EQ.0)THEN
      WRITE(6,*)' Operateur GMV '
      WRITE(6,*)' On ne trouve pas l''indice CENTRE ?'
      GO TO 90
      ENDIF

      CALL LEKTAB(MTAB1,'INCO',INCO)
      IF(INCO.EQ.0)THEN
      WRITE(6,*)' Operateur GMV '
      WRITE(6,*)'Il n''y a pas de table INCO '
      RETURN
      ENDIF
      SEGACT INCO

      CALL KRIPAD(MELEM1,IZIPAD)
      SEGACT MELEMC
      NBC=MELEMC.NUM(/2)
      CALL RSETXI(IZIPAD.LECT,MELEMC.NUM,NBC)
      SEGDES MELEMC
C*****************************************************************************

      CALL LEKTAB(MTAB1,'KIZG',KIZG)
      IF(KIZG.EQ.0)THEN
      CALL CRTABL(KIZG)
      CALL ECMM(KIZG,'SOUSTYPE','KIZG')
      CALL ECMO(MTAB1,'KIZG','TABLE   ',KIZG)
      ELSE
      SEGACT KIZG
      ENDIF

C VERIFICATIONS SUR LES INCONNUES
      NBINC=LINCO.MOTS(/2)
      IF(NBINC.NE.1)THEN
      WRITE(6,*)' Operateur GMV '
      WRITE(6,*)'Nombre d''inconnues incorrecte : ',NBINC,' On attend 1'
      RETURN
      ENDIF

      NOMI=LINCO.MOTS(1)
      SEGACT INCO

      TYPE=' '
      CALL ACMO(INCO,NOMI,TYPE,MCHPOI)
      IF(TYPE.NE.'CHPOINT ')THEN
      WRITE(6,*)' Operateur GMV '
      WRITE(6,*)' L objet CHPOINT ',NOMI,' n existe pas dans la table'
      RETURN
      ELSE
      CALL LICHT(MCHPOI,IZTU1,TYPC,IGEOM0)
      ENDIF

      TYPE=' '
      CALL ACMO(KIZG,NOMI,TYPE,IZG1)
      IF(TYPE.NE.'CHPOINT ')THEN
      NC=IZTU1.VPOCHA(/2)
      TYPE='SOMMET'
      CALL CRCHPT(TYPE,IGEOM0,NC,2,IZG1)
      CALL ECMO(KIZG,NOMI,'CHPOINT ',IZG1)
      ENDIF

      CALL LICHT(IZG1,IZGG1,TYPC,IGEOM)

      TYPE=' '
      CALL ACMO(MTAB1,'PRESSION',TYPE,MTABP)
      IF(TYPE.NE.'TABLE')THEN
      WRITE(6,*)' Operateur GMV '
      WRITE(6,*)' Il n''y a pas de table pression'
      RETURN
      ELSE
      SEGACT MTABP
      TYPE=' '
      CALL ACMO(MTABP,'PRESSION',TYPE,MCHP)
      IF(TYPE.NE.'CHPOINT')THEN
      WRITE(6,*)' Operateur GMV '
      WRITE(6,*)' Il n''y a pas de Champ de pression'
      RETURN
      ELSE
      CALL LICHT(MCHP,IZPP,TYPE,IGEOMP)
      IF(IZPP.EQ.0)GO TO 90
      ENDIF

      ENDIF

      SEGACT MELEME,MENTR
      NE=MENTR.NUM(/2)
      PE=0.
      DO 31 I=1,NE
      I1=MENTR.NUM(1,I)
      NPP1=IZIPAD.LECT(I1)
      PE=PE+IZPP.VPOCHA(NPP1,1)
 31   CONTINUE
      PE=PE/FLOAT(NE)

      SEGACT MSORT
      NS=MSORT.NUM(/2)
      PS=0.
      DO 32 I=1,NS
      I1=MSORT.NUM(1,I)
      NPP1=IZIPAD.LECT(I1)
      PS=PS+IZPP.VPOCHA(NPP1,1)
 32   CONTINUE
      PS=PS/FLOAT(NS)
      DELTAP=PS-PE

      SEGACT MDEBI
      NNP=MDEBI.NUM(/1)
      NNE=MDEBI.NUM(/2)
      DO 33 K=1,NNE
      DO 331 I=1,NNP
      I1=MDEBI.NUM(I,K)
      IF(IZIPAD.LECT(I1).EQ.0)THEN
      WRITE(6,*)' Operateur GMV '
      WRITE(6,*)' LDEBIT SPG incompatible '
      RETURN
      ENDIF
 331  CONTINUE
 33   CONTINUE

      IMPRD=0
      CALL FFDBIT(IZTU1,MDEBI,IZIPAD,IAXI,Q,IMPRD)


      IF(IKIMP.EQ.0)THEN
      SEGACT MQR,MPR
        NPC=MPR.PROG(/1)

        CALL INTEPL(DELTAP,MQR.PROG,MPR.PROG,NPC,QTH,IRET)

      CALL ACMF(MTABA,'K0',AK0)
      IF(QTH.LT.Q)THEN
      AK1=AK0/(1.D0+(Q-QTH)/(Q+1.E-20))
      ELSE
      IF(DELTAP.GE.0)THEN
      AK1=AK0*QTH/(Q+1.E-20)
      ELSE
      AK1=AK0*Q/(QTH+1.E-20)
      ENDIF
      ENDIF
      AK0=W1*AK1+(1.-W1)*AK0

      ELSE
        AK0=AKIMP
      ENDIF

      IF(IMPR.NE.0)THEN
      IMP=MOD(IPAS,IMPR)
      IF(IMP.EQ.0)THEN

      IF(IKIMP.NE.0)THEN
      WRITE(6,*)' GMV : CAS K IMPOSE ',AKIMP
      QTH=0.
      ENDIF
      WRITE(6,1888) DELTAP,Q,QTH,AK0
 1888 FORMAT(1X,' GMV : DELTAP=',1PE12.4,' QR=',1PE11.4,' QTh=',1PE11.4,
     &' AK0=',1PE11.4)
      ENDIF
      ENDIF

      XVEC2(1)=XVEC2(1)*AK0
      XVEC2(2)=XVEC2(2)*AK0
      XVEC2(3)=XVEC2(3)*AK0
      CALL ECMF(MTABA,'K0',AK0)

      NPT=IZGG1.VPOCHA(/1)

      SEGACT MELEME
      NBSOUS=LISOUS(/1)
      IF(NBSOUS.NE.0)THEN
      WRITE(6,*)' Operateur GMV '
      WRITE(6,*)
     & ' La zone doit etre compose d''un seul type d''element'
      RETURN
      ENDIF

      SEGACT IZTU1
      NBSOUS=1
      NUTOEL=0

      DO 1 L=1,NBSOUS
      IPT1=MELEME
      IF(NBSOUS.NE.1)IPT1=LISOUS(L)
      SEGACT IPT1

      NP  =IPT1.NUM(/1)
      NBEL=IPT1.NUM(/2)

      IES=IDIM
C
      CALL ZGMV(IPT1.NUM,NBEL,NUTOEL,NPT,IES,NP,IAXI,
     & IZIPAD.LECT,XVEC2,
     & IZTU1.VPOCHA,
     & IZGG1.VPOCHA,
     & IZVOL.VPOCHA)

      SEGDES IPT1
      NUTOEL=NUTOEL+NBEL

 1    CONTINUE
      SEGDES MELEME

      IF(IKIMP.EQ.0)SEGDES MPR,MQR
      SEGDES IZTU1
      SEGDES IZGG1
      SEGDES IZVOL
      SEGDES LINCO
      SEGDES MTABX,MTAB1,INCO,KIZG,MTABA,MTABP,MTABD,MTABZ
      SEGDES IZPP,MDEBI,MSORT,MENTR
      SEGSUP IZIPAD

 89   CONTINUE
      IPAS=IPAS+1
      RETURN
 90   CONTINUE
      WRITE(6,*)' Interuption anormale de GMV'
      RETURN
 1002 FORMAT(10(1X,1PE11.4))
      END