Télécharger ytscal.eso

Retour à la liste

ytscal
C YTSCAL    SOURCE    FANDEUR   22/01/03    21:16:00     11136          
      SUBROUTINE YTSCAL
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
C***********************************************************************
C
C  CET OPERATEUR DISCRETISE UNE EQUATION DE TRANSPORT
C    EN 2D SUR LES ELEMENTS QUA4 ET TRI3      PLAN OU AXI
C    EN 3D SUR LES ELEMENTS CUB8 ET PRI6
C    LES OPERATEURS SONT "SOUS-INTEGRES"
C
C  SYNTAXE :
C  ---------
C    1/ Cas incompréssible
C
C de/dt + u Grad e = alpha Lapl e  < + S >
C
C   'OPER' 'TSCAL' al 'UN' s          'INCO' EN :
C   'OPER' 'TSCAL' al 'UN' s alt sgt  'INCO' EN :
C
C    2/ Cas compréssible
C
C dh/dt + u Grad h + h Div u = alpha Lapl(tn) < + S >
C
C   'OPER' 'TSCAL' lb 'UN' s tn         'INCO' HN :
C   'OPER' 'TSCAL' lb 'UN' s tn lbt sgt 'INCO' HN :
C
C
C    al,alt Diffusivité thermique moléculaire ou turbulente
C         FLOTTANT où CHPOINT SCAL CENTRE
C    lb,lbt Conductivité thermiquemoléculaire ou turbulente
C         FLOTTANT où CHPOINT SCAL CENTRE
C    sgt  Prandtl turbulent
C    s    source volumique
C         POINT où CHPOINT SCAL CENTRE
C
C    un   Champ de vitesse transportant
C         CHPOINT VECT SOMMET
C    tn   Champ de température
C         CHPOINT SCAL SOMMET
C
C***********************************************************************
 
-INC CCVQUA4
 
-INC PPARAM
-INC CCOPTIO
-INC CCGEOME
-INC SIZFFB
      POINTEUR IZF1.IZFFM,IZH2.IZHR
      PARAMETER (LRV=64)
      SEGMENT PETROV
      REAL*8 WT(LRV,NP,NPG,KDIM),WS(LRV,NP,NPG,KDIM),HK(LRV,IDIM,NP,NPG)
      REAL*8 UIL(LRV,IDIM,NPG),DUIL(LRV,IDIM,NPG)
      REAL*8 PGSK(LRV,NPG),RPGK(LRV,NPG),AIRE(LRV),ANUK(LRV),COEFK(LRV)
      REAL*8 AJK(LRV,IDIM,IDIM,NPG)
      ENDSEGMENT
 
-INC SMCHAML
-INC SMCOORD
-INC SMLENTI
      POINTEUR IPADI.MLENTI,IPADS.MLENTI
      POINTEUR IPADD.MLENTI,IPADQ.MLENTI
-INC SMELEME
      POINTEUR MELEMC.MELEME,MELEMS.MELEME
      POINTEUR MELEMQ.MELEME,MELEMK.MELEME
      POINTEUR MELEMI.MELEME,MELEP1.MELEME
-INC SMCHPOI
      POINTEUR MCHPIN.MCHPOI
      POINTEUR IZGG1.MPOVAL
      POINTEUR IZTU1.MPOVAL
      POINTEUR VISCO.MPOVAL
      POINTEUR IZTGG3.MPOVAL,IZTGG4.MPOVAL
      POINTEUR MZALT.MPOVAL,MZST.MPOVAL
      POINTEUR IZVOL.MPOVAL,IZTCO.MPOVAL
      POINTEUR VITESS.MPOVAL,UTRANS.MPOVAL
-INC SMMATRIK
      POINTEUR IPM.IZAFM
      SEGMENT IMATRS
      INTEGER LIZAFS(NBSOUS,NBME)
      ENDSEGMENT
      POINTEUR IPMS.IZAFM,IPS1.IZAFM,IPS2.IZAFM,IPS3.IZAFM
-INC SMLMOTS
      POINTEUR LINCO.MLMOTS
      CHARACTER*8 NOMZ,NOMI,NOM0,TYPE,TYPC,NOM,NOMA,MTYP,CHAI
      CHARACTER*4 NOM4(3)
      PARAMETER (NTB=1)
      CHARACTER*8 LTAB(NTB)
      DIMENSION KTAB(NTB),IXV(3),RO(1)
      SAVE IPAS
      DATA LTAB/'KIZX    '/,IPAS/0/
C*****************************************************************************
CTSCAL
C      write(6,*)' Debut TSCAL'
 
      segact mcoord
 
C Nouvelle directive EQUA de EQEX
      MTYP=' '
      ro(1)=1.d0
      CALL QUETYP(MTYP,0,IRET)
      IF(IRET.EQ.0)THEN
C% On attend un des objets : %m1:8 %m9:16 %m17:24 %m25:32 %m33:40
      MOTERR( 1: 8) = 'CHAI    '
      MOTERR( 9:16) = 'MMODEL  '
      MOTERR(17:24) = 'TABLE   '
      CALL ERREUR(472)
      RETURN
      ENDIF
 
      IF(MTYP.EQ.'MMODEL')THEN
        CALL YTCLSF('TCLS    ','TSCA    ')
        RETURN
      ELSEIF(MTYP.EQ.'MOT     ')THEN
        CALL LIRCHA(CHAI,1,IRET)
        CALL YTCLSF(CHAI,'TSCA    ')
        RETURN
      ENDIF
C Fin Nouvelle directive EQUA de EQEX
 
c     CALL LITABS(LTAB,KTAB,NTB,1,IRET)
c     IF(IERR.NE.0) RETURN
c     MTABX=KTAB(1)
 
C
C- Récupération de la table xNS type KIZX (pointeur MTABX)
C
C     write(6,*)'Récupération de la table xOP type KIZX pointeur MTABX'
      CALL LIROBJ('TABLE',MTABX,1,IRET)
      IF(IRET.EQ.0)THEN
C% On ne trouve pas d'objet de type %m1:8
      MOTERR( 1: 8) = 'TABLE'
      CALL ERREUR(37)
      RETURN
      ENDIF
C A tout hazard on regarde si une table ne peut en cacher une autre
C     write(6,*)' A tout hazard on regarde si une table ne peut en',
C    &' cacher une autre !!!!!'
 
      MTYP='TABLE'
      CALL QUETYP(MTYP,0,IRET)
      IF(IRET.EQ.1)THEN
      CALL LIROBJ('TABLE',MTABX,1,IRET)
       IF(IRET.EQ.0)THEN
C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24
         MOTERR( 1: 8) = '  MTABX '
         MOTERR( 9:16) = '  MTABX '
         MOTERR(17:24) = '  KIZX  '
         CALL ERREUR(786)
         RETURN
       ENDIF
C     write(6,*)' EH BEN C EST LE CAS une KIZX // !!!!!!!!!!'
cc      call pplist(MTABX)
      TYPE=' '
      CALL ACMO(MTABX,'DOMZ',TYPE,MMDZ)
      CALL LEKMOD(MMDZ,MTABZ,INEFMD)
      TYPE='TABLE'
      CALL ECMO(MTABX,'TDOMZ',TYPE,MTABZ)
      ELSE
      TYPE=' '
      CALL ACMO(MTABX,'DOMZ',TYPE,MMDZ)
       IF(TYPE.EQ.'MMODEL')THEN
c     call pplist(mtabx)
       CALL LEKMOD(MMDZ,MTABZ,INEFMD)
       ENDIF
      ENDIF
 
C.......................................................................
C
C- Récupération de la table RV type EQEX (pointeur MTAB1)
C
      CALL LEKTAB(MTABX,'EQEX',MTAB1)
      IF(MTAB1.EQ.0)THEN
C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24
         MOTERR( 1: 8) = '  EQEX  '
         MOTERR( 9:16) = '  EQEX  '
         MOTERR(17:24) = '  KIZX  '
         CALL ERREUR(786)
         RETURN
      ENDIF
      CALL ACME(MTAB1,'NAVISTOK',NASTOK)
      IF(NASTOK.EQ.0)THEN
      CALL ZTSCAL(MTABX,MTAB1)
      RETURN
      ENDIF
C
C- Récupération de la table INCO (pointeur KINC)
C
      CALL LEKTAB(MTAB1,'INCO',KINC)
      IF(KINC.EQ.0)THEN
C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24
         MOTERR( 1: 8) = '  INCO  '
         MOTERR( 9:16) = '  INCO  '
         MOTERR(17:24) = '  EQEX  '
         CALL ERREUR(786)
      RETURN
      ENDIF
C.......................................................................
 
      CALL ACMM(MTABX,'NOMZONE',NOMZ)
 
C On récupère le MODEL partitionné ou non
      TYPE=' '
      CALL ACMO(MTABX,'DOMZ',TYPE,MMDZ)
       IF(TYPE.NE.'MMODEL')THEN
C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24
         MOTERR( 1: 8) = '  DOMZ  '
         MOTERR( 9:16) = '  DOMZ  '
         MOTERR(17:24) = '  KIZX  '
         CALL ERREUR(786)
         RETURN
       ENDIF
 
C*****************************************************************************
C OPTIONS
C KIMPL = 0 -> EXPL     1 -> IMPL     2 -> SEMI
C KFORM = 0 -> SI       1 -> EF       2 -> VF    3 -> EFMC
C IDCEN = 0-> rien  1-> CENTREE  2-> SUPGDC  3-> SUPG  4-> TVISQUEU  5-> CNG
 
      IAXI=0
      IF(IFOMOD.EQ.0)IAXI=2
C
C- Récupération de la table des options KOPT (pointeur KOPTI)
C
      CALL LEKTAB(MTABX,'KOPT',KOPTI)
      IF (KOPTI.EQ.0) THEN
C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24
         MOTERR( 1: 8) = '  KOPT  '
         MOTERR( 9:16) = '  KOPT  '
         MOTERR(17:24) = '  KIZX  '
         CALL ERREUR(786)
         RETURN
      ENDIF
 
 
      CALL ACME(KOPTI,'MTRMASS ',MMPG)
      IPG=0
      IF(MMPG.EQ.3)IPG=1
      KDIM=1
      CALL ACME(KOPTI,'IDCEN',IDCEN)
      IF(IDCEN.EQ.2)KDIM=IDIM
      CALL ACME(KOPTI,'IKOMP',IKOMP)
      CALL ACME(KOPTI,'KIMPL',KIMPL)
      CALL ACME(KOPTI,'IDIV',IDIV)
      CALL ACMF(KOPTI,'CMD',CMD)
      CALL ACMF(KOPTI,'AIMPL',AIMPL)
      AG=AIMPL
      AD=AIMPL-1.D0
      CALL ACME(KOPTI,'KFORM',KFORM)
      CALL ACME(KOPTI,'KPOIN',KPRE)
      CALL ACME(KOPTI,'KMACO',KMACO)
         IF (IERR.NE.0) RETURN
 
C*****************************************************************************
C
C- Récupération de la table DOMAINE associée au domaine local
C
 
C E/  MMODEL : Pointeur de la table contenant l'information cherchée
C  /S IPOINT : Pointeur sur la table DOMAINE
C  /S INEFMD : Type formulation INEFMD=1 LINE,=2 MACRO,=3 QUADRATIQUE
C                               INEFMD=4 LINB
 
c??      CALL LEKMOD(MMDZ,MTABZ,INEFMD)
 
      CALL LEKTAB(MTABZ,'MAILLAGE',MELEME)
      CALL LEKTAB(MTABZ,'SOMMET',MELEMS)
      CALL LEKTAB(MTABZ,'CENTRE',MELEMC)
      CALL LEKTAB(MTABZ,'MACRO',MACRO)
      CALL LEKTAB(MTABZ,'MACRO1',MELEMK)
      CALL LEKTAB(MTABZ,'QUADRATI',MQUAD)
         IF (IERR.NE.0) RETURN
 
      MELEMQ=MELEMC
      MELEP1=MELEMC
      IF(KPRE.EQ.3)THEN
      CALL LEKTAB(MTABZ,'CENTREP0',MELEMQ)
      MELEP1=MELEMQ
      ELSEIF(KPRE.EQ.4)THEN
      CALL LEKTAB(MTABZ,'CENTREP1',MELEMQ)
      CALL LEKTAB(MTABZ,'ELTP1NC ',MELEP1)
      ENDIF
 
         IF (IERR.NE.0) RETURN
C***
 
 
C*************************************************************************
C VERIFICATIONS SUR LES INCONNUES
C
C- Récupération du nombre d'inconnues et du nom de l'inconnue NOMI
C
      TYPE='LISTMOTS'
      CALL ACMO(MTABX,'LISTINCO',TYPE,LINCO)
      IF (IERR.NE.0) RETURN
      SEGACT LINCO
      NBINC=LINCO.MOTS(/2)
      IF(NBINC.NE.1)THEN
C        Indice %m1:8 : contient plus de %i1 %m9:16
         MOTERR( 1:8) = 'LISTINCO'
         INTERR(1) = 1
         MOTERR(9:16) = ' MOTS   '
         CALL ERREUR(799)
         RETURN
      ENDIF
 
      NOMI=LINCO.MOTS(1)
      NOMA=NOMI
      WRITE(NOM4(1),FMT='(A4)')NOMI(1:4)
C
C- Récupération de l'inconnue
C
      TYPE=' '
      CALL ACMO(KINC,NOMI,TYPE,MCHPOI)
      IF(TYPE.NE.'CHPOINT ')THEN
C        Indice %m1:8 : ne contient pas un objet de type %m9:16
         MOTERR( 1: 8) = 'INC '//NOMI
         MOTERR( 9:16) = 'CHPOINT '
         CALL ERREUR(800)
         RETURN
      ELSE
      CALL LICHTL(MCHPOI,IZTU1,TYPC,MELEMI)
         MCHPIN=MCHPOI
         NINKO = IZTU1.VPOCHA(/2)
         NPTI  = IZTU1.VPOCHA(/1)
         IF (NINKO.NE.1) THEN
C        Indice %m1:8 : Le %m9:16 n'a pas le bon nombre de composantes
            MOTERR( 1: 8) = 'INC '//NOMI
            MOTERR( 9:16) = 'CHPOINT '
            CALL ERREUR(784)
            RETURN
         ENDIF
      ENDIF
 
C*****************************************************************************
C Le domaine de definition est donne par le SPG  de la premiere inconnue
C Les inconnues suivantes devront posseder ce meme pointeur
C On verifie que les points de la zone sont tous inclus dans ce SPG
 
      CALL KRIPAD(MELEMI,IPADI)
      IPADS=IPADI
      NPTS=NPTI
      IPADD=IPADI
      IF(MELEMI.NE.MELEMS)THEN
      CALL KRIPAD(MELEMS,IPADS)
      NPTS=MELEMS.NUM(/2)
      ENDIF
 
      IF(IPAS.EQ.0)THEN
      CALL VERPAD(IPADI,MELEME,IRET)
      IF(IRET.NE.0)THEN
C     Indice %m1:8 : L'objet %m9:16 n'a pas le bon support géométrique
         MOTERR(1: 8) = 'INC '//NOMI
         MOTERR(9:16) = 'CHPOINT '
         CALL ERREUR(788)
         RETURN
      ENDIF
      ENDIF
 
C*****************************************************************************
      IF(KFORM.EQ.1)THEN
      IHV=0
 
        CMT=CMD
        DT=0.D0
        IF(IDCEN.EQ.4.OR.IDCEN.EQ.5)THEN
C IDCEN 1-> CENTREE  2-> SUPGDC  3-> SUPG  4-> TVISQUEU  5-> CNG
        TYPE=' '
        CALL ACMO(MTABT,'DELTAT',TYPE,IENT)
          IF(TYPE.NE.'ENTIER')THEN
            WRITE(6,*)' Opérateur TSCA'
            WRITE(6,*)' On reclame un pas de temps'
C        Option %m1:8 incompatible avec les données
            MOTERR( 1: 8) = '  EF    '
            CALL ERREUR(803)
            RETURN
          ELSE
            CALL ACMF(MTABT,'DELTAT',DT)
            CMT=DT
          ENDIF
        ENDIF
 
      CALL YCLS('TSCA    ',MTABX,MTABZ,IHV,MCHPIN,NOMI,MATRIK,MCHPO1)
 
      CALL ECROBJ('MATRIK',MATRIK)
      CALL ECROBJ('CHPOINT',MCHPO1)
 
      RETURN
      ENDIF
 
C*****************************************************************************
C Lecture des coefficient
C Type du coefficient :
C IK1=0 CHPOINT IK1=1 scalaire  IK1=2 vecteur
 
C     write(6,*)' Opérateur TSCAL lecture des coefficients'
      CALL ACME(MTABX,'IARG',IARG)
      IF(IKOMP.EQ.0)THEN
      IF(IARG.NE.3.AND.IARG.NE.4.AND.IARG.NE.5)THEN
      WRITE(6,*)' Opérateur TSCAL : option non conservative '
      WRITE(6,*)' Nombre d''arguments incorrect : ',IARG
      WRITE(6,*)' On attend 3 , 4 ou 5 '
C           Indice %m1:8 : nombre d'arguments incorrect
            MOTERR(1:8) = 'IARG    '
            CALL ERREUR(804)
            RETURN
         ENDIF
      ELSEIF(IKOMP.EQ.1.AND.KFORM.NE.1)THEN
      IF(IARG.NE.4.AND.IARG.NE.6)THEN
      WRITE(6,*)' Opérateur TSCAL : option compréssible '
      WRITE(6,*)' Nombre d''arguments incorrect : ',IARG
      WRITE(6,*)' On attend 4 ou 6 '
C           Indice %m1:8 : nombre d'arguments incorrect
            MOTERR(1:8) = 'IARG    '
            CALL ERREUR(804)
            RETURN
         ENDIF
      ENDIF
 
 
C--Cas incompressible
      IF(IKOMP.EQ.0.OR.IKOMP.EQ.2)THEN
C 1er coefficient  Alpha
      IXV(1)=MELEMC
      IXV(2)=1
      IXV(3)=0
      CALL LEKCOF('Opérateur TSCAL :',
     & MTABX,KINC,1,IXV,MCOF,VISCO,NPT1,NC1,IKL,IRET)
      IF(IRET.EQ.0)RETURN
 
      IK2=-1
      IKG=-1
      IK3=-1
      IK4=-1
      IZTGG4=IZTU1
      MZALT=VISCO
      MZST=VISCO
 
C 2ème coefficient  UN
      IXV(1)=-MELEMS
      IXV(2)=0
      IXV(3)=1
      CALL LEKCOF('Opérateur TSCAL :',
     & MTABX,KINC,2,IXV,MUTR,UTRANS,NPTU,NC2,IKU,IRET)
      IF(IRET.EQ.0)RETURN
      IF(IKU.EQ.2)IKU=1
 
C 3ème coefficient  source
      IXV(1)=MELEMC
      IXV(2)=1
      IXV(3)=0
      CALL LEKCOF('Opérateur TSCAL :',
     & MTABX,KINC,3,IXV,IZTG3,IZTGG3,NPT3,NC3,IKS,IRET)
      IF(IRET.EQ.0)RETURN
 
      IF(IARG.EQ.5)THEN
C alpha turbulent
      IXV(1)=MELEMC
      IXV(2)=0
      IXV(3)=0
      CALL LEKCOF('Opérateur TSCAL :',
     & MTABX,KINC,4,IXV,MALT,MZALT,NPT4,NC4,IK4,IRET)
      IF(IRET.EQ.0)RETURN
 
C sigma turbulent
      IXV(1)=0
      IXV(2)=1
      IXV(3)=0
      CALL LEKCOF('Opérateur TSCAL :',
     & MTABX,KINC,5,IXV,MST,MZST,NPT5,NC5,IK5,IRET)
      IF(IRET.EQ.0)RETURN
      IF(MZST.VPOCHA(1,1).LE.0.D0)THEN
      WRITE(6,*)' Opérateur TSCAL :'
      WRITE(6,*)'Valeur du Prandtl turbulent érronée'
      RETURN
      ENDIF
      ENDIF
 
C--Cas compréssible
      ELSEIF(IKOMP.EQ.1)THEN
C 1er coefficient  Lambda
      IXV(1)=MELEMC
      IXV(2)=1
      IXV(3)=0
      CALL LEKCOF('Opérateur TSCAL :',
     & MTABX,KINC,1,IXV,MCOF,VISCO,NPT1,NC1,IKL,IRET)
      IF(IRET.EQ.0)RETURN
 
      IK2=-1
      IKG=-1
      IK3=-1
      IK4=-1
      MZALT=VISCO
      MZST=VISCO
 
C 2ème coefficient  UN
      IXV(1)=-MELEMS
      IXV(2)=0
      IXV(3)=1
      CALL LEKCOF('Opérateur TSCAL :',
     & MTABX,KINC,2,IXV,MUTR,UTRANS,NPTU,NC2,IKU,IRET)
      IF(IRET.EQ.0)RETURN
      IF(IKU.EQ.2)IKU=1
 
C 3ème coefficient  source
      IXV(1)=MELEMC
      IXV(2)=1
      IXV(3)=0
      CALL LEKCOF('Opérateur TSCAL :',
     & MTABX,KINC,3,IXV,IZTG3,IZTGG3,NPT3,NC3,IKS,IRET)
      IF(IRET.EQ.0)RETURN
 
C 4ème coefficient  tn
      IXV(1)=MELEMS
      IXV(2)=0
      IXV(3)=0
      CALL LEKCOF('Opérateur TSCAL :',
     & MTABX,KINC,4,IXV,IZTG4,IZTGG4,NPT4,NC4,IK4,IRET)
      IPADD=IPADS
      IF(IRET.EQ.0)RETURN
 
      IF(IARG.EQ.6)THEN
C lambda turbulent
      IXV(1)=MELEMC
      IXV(2)=0
      IXV(3)=0
      CALL LEKCOF('Opérateur TSCAL :',
     & MTABX,KINC,5,IXV,MALT,MZALT,NPT5,NC5,IK5,IRET)
      IF(IRET.EQ.0)RETURN
 
C sigma turbulent
      IXV(1)=0
      IXV(2)=1
      IXV(3)=0
      CALL LEKCOF('Opérateur TSCAL :',
     & MTABX,KINC,6,IXV,MST,MZST,NPT6,NC6,IK6,IRET)
      IF(IRET.EQ.0)RETURN
      IF(MZST.VPOCHA(1,1).LE.0.D0)THEN
        call erreur(992)
      RETURN
      ENDIF
      ENDIF
 
      ENDIF
 
C     write(6,*)' Opérateur TSCAL : Fin lecture Arguments '
C Fin lecture Arguments ************************************************
 
C*****************************************************************************
C Lecture table 'PASDETPS'
      CALL LEKTAB(MTAB1,'PASDETPS',MTABT)
      IF(MTABT.EQ.0)THEN
      CALL CRTABL(MTABT)
      CALL ECMM(MTABT,'SOUSTYPE','PASDETPS')
      CALL ECMO(MTAB1,'PASDETPS','TABLE   ',MTABT)
      DT=1.D30
      DTP=1.D30+DT
      IPAT=1
      CALL ECME(MTABT,'NUPASDT',IPAT)
      DTM1=1.D-20
      CALL ECMF(MTABT,'DELTAT-1',DTM1)
      ELSE
      CALL ACMF(MTABT,'DELTAT',DTP)
      CALL ACMF(MTABT,'DELTAT-1',DTM1)
      ENDIF
 
 
C*********** FORMULATIONS **********
C*********** FORMULATIONS **********
C*********** FORMULATIONS **********
 
 
 
C************************** FORMULATION EFM1 ***************************
      IF(KFORM.EQ.0)THEN
C Formulation EFM1
C Vérification des options
 
      IF(INEFMD.NE.1.AND.INEFMD.NE.2)THEN
C        Option %m1:8 incompatible avec les données
            MOTERR( 1: 8) = '  EFM1  '
      CALL ERREUR(803)
      RETURN
      ENDIF
 
      IF(KIMPL.NE.0)THEN
C        Option %m1:8 incompatible avec les données
            MOTERR( 1: 8) = ' EFM1   '
            WRITE(6,*)' Options incompatibles : EFM1 et (IMPL ou SEMI) '
            CALL ERREUR(803)
            RETURN
        ENDIF
 
      CALL LEKTAB(MTABZ,'XXPSOML',MCHELM)
      CALL LEKTAB(MTABZ,'XXDXDY',MCHPDX)
      CALL LEKTAB(MTABZ,'XXVOLUM',MCHPVO)
         IF (IERR.NE.0) RETURN
 
      SEGACT MCHELM
      CALL LICHTL(MCHPDX,IZTCO,TYPC,IGEOM)
      NELZ=IZTCO.VPOCHA(/1)
      CALL LICHTL(MCHPVO,IZVOL,TYPC,IGEOM)
 
      IF(IKOMP.EQ.1.AND.(IDCEN.EQ.6.OR.IDCEN.EQ.7))THEN
C        Option %m1:8 incompatible avec les données
            MOTERR( 1: 8) = ' EFM1   '
            WRITE(6,*)' Option de décentrement non prévue en',
     &                ' formulation conservative '
            CALL ERREUR(803)
            RETURN
        ENDIF
 
      IF(IDCEN.GE.4.AND.(IDCEN.NE.6.AND.IDCEN.NE.7))THEN
C        Option %m1:8 incompatible avec les données
            MOTERR( 1: 8) = ' EFM1   '
            WRITE(6,*)' Option de décentrement non prévue en',
     &                ' formulation EFM1 '
            CALL ERREUR(803)
            RETURN
        ENDIF
 
      IF(IDIM.EQ.3.AND.(IDCEN.EQ.6.OR.IDCEN.EQ.7))THEN
C        Option %m1:8 incompatible avec les données
            MOTERR( 1: 8) = ' EFM1   '
            WRITE(6,*)' Option de décentrement non prévue en 3D'
            CALL ERREUR(803)
            RETURN
        ENDIF
 
      NC=IZTU1.VPOCHA(/2)
      TYPE='SOMMET'
      CALL KRCHPT(TYPE,MELEMS,NC,IZG1,NOM4)
 
      CALL LICHTM(IZG1,IZGG1,TYPC,IGEOM)
 
      NPTS=IZGG1.VPOCHA(/1)
 
 
      CALL LEKTAB(MTABZ,'MATESI',MATRIK)
      IF (IERR.NE.0) RETURN
      SEGACT MATRIK
 
      IMATRI=IRIGEL(4,1)
      SEGACT IMATRI
C---
 
      SEGACT MELEME
      NBSOUS=LISOUS(/1)
      IF(NBSOUS.EQ.0)NBSOUS=1
      NUTOEL=0
      DT=1.E30
 
      DO 1 L=1,NBSOUS
      IPT1=MELEME
      IF(NBSOUS.NE.1)IPT1=LISOUS(L)
      SEGACT IPT1
 
      NOM0=NOMS(IPT1.ITYPEL)//'    '
      CALL KALPBG(NOM0,'FONFORM0',IZFFM)
      SEGACT IZFFM*MOD
 
      MCHAML=ICHAML(L)
      SEGACT MCHAML
      MELVAL=IELVAL(1)
      SEGACT MELVAL
 
      IF(IMACHE(L).NE.IPT1)THEN
       write(*,*)'IPT1,IMACHE ',IPT1,IMACHE(L)
       goto 90
      ENDIF
 
      IZAFM=LIZAFM(L,1)
      IPM1=IZAFM
      SEGACT IZAFM
      IF(IAXI.NE.0)THEN
      IPM1=LIZAFM(L,2)
      SEGACT IPM1
      ENDIF
 
      NP  =IPT1.NUM(/1)
      NBEL=IPT1.NUM(/2)
      IES=IDIM
 
      NPTT=IZTGG4.VPOCHA(/1)
      NPTI=IZTGG4.VPOCHA(/1)
 
      IF(IDCEN.GE.1.AND.IDCEN.LE.3)THEN
 
      IF(IKOMP.EQ.0)THEN
 
      CALL YCTSCL(AM,IPM1.AM,VELCHE,IPT1.NUM,NBEL,NUTOEL,IES,NP,
     & IAXI,IPADI.LECT,IPADS.LECT,IPADD.LECT,IKOMP,IARG,
     & VISCO.VPOCHA,IKL,UTRANS.VPOCHA,IKU,NPTS,IZTGG4.VPOCHA,NPTT,
     & IZTGG3.VPOCHA,IKS,
     & IZTU1.VPOCHA,IZGG1.VPOCHA,NPTI,
     & MZALT.VPOCHA,MZST.VPOCHA,
     & IZVOL.VPOCHA,IZTCO.VPOCHA,NELZ,IDCEN,IPG,
     & DTM1,DT,DTT1,DTT2,NUEL,DIAEL,IZFFM.FN)
 
      ELSEIF(IKOMP.EQ.1)THEN
 
      CALL YCTSCA(AM,IPM1.AM,VELCHE,IPT1.NUM,NBEL,NUTOEL,NPT,IES,NP,
     & IAXI,IPADI.LECT,IKOMP,IARG,IPADS.LECT,IPADD.LECT,
     & VISCO.VPOCHA,IKL,UTRANS.VPOCHA,IKU,NPTS,IZTGG4.VPOCHA,
     & IZTGG3.VPOCHA,IKS,IZTU1.VPOCHA,
     & IZGG1.VPOCHA,MZALT.VPOCHA,MZST.VPOCHA,
     & IZVOL.VPOCHA,IZTCO.VPOCHA,NELZ,IDCEN,
     & DTM1,DT,DTT1,DTT2,NUEL,DIAEL)
 
      ENDIF
 
      ELSEIF(IDCEN.EQ.6)THEN
 
 
      N=NPTS
      NC=1
      SEGINI MPOVA6
 
      CALL YPSI(AM,IPM1.AM,VELCHE,IPT1.NUM,NBEL,NUTOEL,NPT,IES,NP,
     & IAXI,IPADI.LECT,IKOMP,IARG,
     & VISCO.VPOCHA,IKL,UTRANS.VPOCHA,IKU,NPTS,IPADS.LECT,
     & IZTGG4.VPOCHA,IZTGG3.VPOCHA,IKS,IZTU1.VPOCHA,
     & IZGG1.VPOCHA,MZALT.VPOCHA,MZST.VPOCHA,
     & IZVOL.VPOCHA,IZTCO.VPOCHA,NELZ,MPOVA6.VPOCHA,
     & DTM1,DT,DTT1,DTT2,NUEL,DIAEL)
 
      SEGSUP MPOVA6
 
 
 
      ELSEIF(IDCEN.EQ.7)THEN
 
      CALL YJOHNS(AM,IPM1.AM,VELCHE,IPT1.NUM,NBEL,NUTOEL,NPT,IES,NP,
     & IAXI,IPADI.LECT,IKOMP,IARG,
     & VISCO.VPOCHA,IKL,UTRANS.VPOCHA,IKU,NPTS,IPADS.LECT,
     & IZTGG4.VPOCHA,IZTGG3.VPOCHA,IKS,IZTU1.VPOCHA,
     & IZGG1.VPOCHA,MZALT.VPOCHA,MZST.VPOCHA,
     & IZVOL.VPOCHA,IZTCO.VPOCHA,NELZ,
     & DTM1,DT,DTT1,DTT2,NUEL,DIAEL)
 
      ENDIF
 
 
      SEGDES IZAFM*NOMOD,IPT1*NOMOD,MCHAML*NOMOD,MELVAL*NOMOD
      IF(IAXI.NE.0)SEGDES IPM1*NOMOD
      NUTOEL=NUTOEL+NBEL
 
 1    CONTINUE
      SEGDES MATRIK*NOMOD,MCHELM*NOMOD
      SEGDES IMATRI
 
      IF(DT.LT.DTP)THEN
      CALL ECMF(MTABT,'DELTAT',DT)
      CALL ECMM(MTABT,'OPER','TSCAL')
      CALL ECMM(MTABT,'ZONE',NOMZ)
      CALL ECMF(MTABT,'DTCONV',DTT1)
      CALL ECMF(MTABT,'DTDIFU',DTT2)
      CALL ECMF(MTABT,'DIAEL',DIAEL)
      CALL ECME(MTABT,'NUEL',NUEL)
      ENDIF
      SEGDES VISCO,UTRANS
      IF(IKOMP.EQ.0.AND.IARG.EQ.3)SEGDES IZTGG3*NOMOD
      IF(IKOMP.EQ.1.AND.IARG.EQ.4)SEGDES IZTGG3*NOMOD,IZTGG4*NOMOD
      IF(IKOMP.EQ.1.AND.IARG.EQ.6)SEGDES IZTGG3*NOMOD,IZTGG4*NOMOD
     & ,MZALT*NOMOD,MZST*NOMOD
      SEGDES IZTU1
      SEGDES IZGG1
      SEGDES IZVOL,IZTCO
      SEGDES LINCO
 
      NRIGE=7
      NKID =9
      NKMT =7
      NMATRI=0
      SEGINI MATRIK
      CALL ECROBJ('MATRIK',MATRIK)
      CALL ECROBJ('CHPOINT',IZG1)
 
C*************************************************************************
      ELSE IF(KFORM.EQ.1)THEN
C CAS FORMULATION EF
 
      DT=0.D0
      IF(IDCEN.EQ.4)THEN
      TYPE=' '
      CALL ACMO(MTABT,'DELTAT',TYPE,IENT)
      IF(TYPE.NE.'ENTIER')THEN
      WRITE(6,*)' Opérateur TSCA'
      WRITE(6,*)' On reclame un pas de temps'
C        Option %m1:8 incompatible avec les données
            MOTERR( 1: 8) = '  EF    '
            CALL ERREUR(803)
      RETURN
      ELSE
      CALL ACMF(MTABT,'DELTAT',DT)
      ENDIF
      ENDIF
 
      IHV=0
      NUTOEL=0
      NINKO=IZTU1.VPOCHA(/2)
 
      TYPE=' '
      CALL ACMO(MTABX,'MATELM',TYPE,MATRIK)
      IF(TYPE.EQ.'MATRIK'.AND.KMACO.NE.0)THEN
      SEGACT MATRIK
      NMATRI=IRIGEL(/2)
      MELEME=IRIGEL(1,1)
      SEGACT MELEME
      IMATRI=IRIGEL(4,1)
      SEGACT IMATRI
      NBME=LIZAFM(/2)
      NINKO=NBME
      MELEMS=KSPGP
 
 
      ELSE
 
      NRIGE=7
      NKID =9
      NKMT =7
      NMATRI=1
      SEGINI MATRIK
      IRIGEL(1,1)=MELEME
      IRIGEL(2,1)=MELEME
      IRIGEL(7,1)=2
      NBME=NINKO
      NBELC=0
      CALL LEKTAB(MTABZ,'SOMMET',MELEMS)
 
C --- CALCUL DU NOMBRE DE PAQUETS DE LRV ELEMENTS ---
      CALL KPLRVM(MELEME,LRV,NBSOUS)
 
      SEGINI IMATRI,IMATRS
      IRIGEL(4,1)=IMATRI
      KSPGP=MELEMS
      KSPGD=MELEMS
 
      IF(NBME.EQ.1)THEN
      LISPRI(1)=NOMI(1:4)//'    '
      LISDUA(1)=NOMA(1:4)//'    '
      ELSE
      DO 102 I=1,NBME
      WRITE(NOM,FMT='(I1,A7)')I,NOMI(1:7)
      LISPRI(I)=NOM(1:4)//'    '
      WRITE(NOM,FMT='(I1,A7)')I,NOMA(1:7)
      LISDUA(I)=NOM(1:4)//'    '
 102  CONTINUE
      ENDIF
 
      SEGACT MELEME
      NBSOU1=LISOUS(/1)
      IF(NBSOU1.EQ.0)NBSOU1=1
      LL=0
      NUTOEL=0
      DO 101 L=1,NBSOU1
      IPT1=MELEME
      IF(NBSOU1.NE.1)IPT1=LISOUS(L)
      SEGACT IPT1
      NOM0=NOMS(IPT1.ITYPEL)//'    '
      CALL KALPBG(NOM0,'FONFORM ',IZFFM)
      SEGACT IZFFM*MOD
      IZHR=KZHR(1)
      IZH2=KZHR(2)
      SEGACT IZHR*MOD
      NES=GR(/1)
      NPG=GR(/3)
      IZF1=KTP(1)
 
      NP = IPT1.NUM(/1)
      MP = NP
      NBELEM=IPT1.NUM(/2)
 
      NNN=MOD(NBELEM,LRV)
      IF(NNN.EQ.0) NPACK=NBELEM/LRV
      IF(NNN.NE.0) NPACK=1+(NBELEM-NNN)/LRV
 
      DO 701 KPACK=1,NPACK
 
C --- POUR CHAQUE PAQUET DE LRV ELEMENTS ou moins
         KDEB=1+(KPACK-1)*LRV
         KFIN=MIN(NBELEM,KDEB+LRV-1)
 
      NBEL=KFIN-KDEB+1
      LL=LL+1
 
      SEGINI IPM1,IPS1
      LIZAFM(LL,1)=IPM1
      LIZAFS(LL,1)=IPS1
      IPM2=IPM1
      IPM3=IPM1
      IPS2=IPS1
      IPS3=IPS1
      IF(NBME.GE.2)THEN
      SEGINI IPM2,IPS2
      LIZAFM(LL,2)=IPM2
      LIZAFS(LL,2)=IPS2
      ENDIF
      IF(NBME.GE.3)THEN
      SEGINI IPM3,IPS3
      LIZAFM(LL,3)=IPM3
      LIZAFS(LL,3)=IPS3
      ENDIF
 
 
      KITT=2
      KJTT=IKL
 
      SEGINI PETROV
 
      CALL XCONV(FN,GR,PG,XYZ,HR,PGSQ,RPG,AJ,
     & NES,IDIM,NP,NPG,IAXI,AG,AD,IDIV,CMD,IKOMP,LRV,
     & WT,WS,HK,PGSK,RPGK,AIRE,AJK,UIL,DUIL,COEFK,ANUK,
     & RO,1,UTRANS.VPOCHA,NPTS,IPADS.LECT,VISCO.VPOCHA,IKL,
     & IPT1.NUM(1,KDEB),NBEL,NUTOEL,XCOOR,
     & IPM1.AM,IPM2.AM,IPM3.AM,
     & IPS1.AM,IPS2.AM,IPS3.AM,
     & NINKO,IDCEN,DT,
     & IZTU1.VPOCHA,NPTG,IPADI.LECT)
 
      SEGSUP PETROV
 
      CALL XLAPL(FN,GR,PG,XYZ,HR,PGSQ,RPG,NES,IDIM,NP,NPG,IAXI,
     & VISCO.VPOCHA,VISCO.VPOCHA,VISCO.VPOCHA,KITT,KJTT,IKL,
     & IPT1.NUM(1,KDEB),NBEL,NUTOEL,XCOOR,AIMPL,IKOMP,
     & IPM1.AM,IPM2.AM,IPM3.AM,
     & IPS1.AM,IPS2.AM,IPS3.AM,
     & NINKO,IHV,IARG,VISCO.VPOCHA)
 
      NUTOEL=NUTOEL+NBEL
      SEGDES IPM1,IPS1
      IF(NBME.GE.2)SEGDES IPM2,IPS2
      IF(NBME.GE.3)SEGDES IPM3,IPS3
 
 701  CONTINUE
      SEGDES IPT1
      SEGSUP IZFFM,IZHR,IZF1,IZH2
 101  CONTINUE
      IF(NBSOU1.NE.1)SEGDES MELEME
      ENDIF
 
C On cree le CHPOINT pour les sources et eventuellement le semi
C
      NAT=2
      NSOUPO=1
      SEGACT MELEMS
      N=MELEMS.NUM(/2)
      NC=NINKO
      SEGINI MCHPO1,MSOUP1,MPOVA1
      MCHPO1.IFOPOI=IFOUR
      MCHPO1.MOCHDE=TITREE
      MCHPO1.MTYPOI='SMBR'
      MCHPO1.JATTRI(1)=2
      MCHPO1.IPCHP(1)=MSOUP1
      DO 177 N=1,NINKO
      MSOUP1.NOCOMP(N)=LISDUA(N)
 177  CONTINUE
      MSOUP1.IGEOC=MELEMS
      MSOUP1.IPOVAL=MPOVA1
      SEGDES MCHPO1,MSOUP1
 
      CALL KRIPAD(MELEMQ,IPADQ)
      SEGDES MELEMQ,MELEMS
 
      SEGACT MELEMK
      NBSOUS=MELEMK.LISOUS(/1)
      IF(NBSOUS.EQ.0)NBSOUS=1
 
      NUTOEL=0
      DO 1102 L=1,NBSOUS
      IPT1=MELEMK
      IF(NBSOUS.NE.1)IPT1=MELEMK.LISOUS(L)
      SEGACT IPT1
 
       NOM0=NOMS(IPT1.ITYPEL)//'    '
C      write(6,*)'MQUAD,MACRO,INEFMD=',MQUAD,MACRO,INEFMD
       IF(INEFMD.EQ.3)THEN
C QUAD
       IF(KPRE.EQ.2)NOM0=NOMS(IPT1.ITYPEL)//'PRP0'
       IF(KPRE.EQ.3)NOM0=NOMS(IPT1.ITYPEL)//'PRP0'
       IF(KPRE.EQ.4)NOM0=NOMS(IPT1.ITYPEL)//'PRP1'
       ELSEIF(INEFMD.EQ.2)THEN
C MACRO
       IF(KPRE.EQ.2)NOM0=NOMS(IPT1.ITYPEL)//'    '
       IF(KPRE.EQ.3)NOM0=NOMS(IPT1.ITYPEL)//'MCP0'
       IF(KPRE.EQ.4)NOM0=NOMS(IPT1.ITYPEL)//'MCP1'
       ELSE
       IF(KPRE.EQ.2)NOM0=NOMS(IPT1.ITYPEL)//'    '
       ENDIF
 
       CALL KALPBG(NOM0,'FONFORM ',IZFFM)
 
 
      SEGACT IZFFM*MOD
      IZHR=KZHR(1)
      IZH2=KZHR(2)
      SEGACT IZHR*MOD
      NES=GR(/1)
      NPG=GR(/3)
      IZF1=KTP(1)
      SEGACT IZF1*MOD
      MP1=IZF1.FN(/1)
      NP = IPT1.NUM(/1)
      MP = NP
      NBEL=IPT1.NUM(/2)
      NELG=IZTGG3.VPOCHA(/1)
      NPT =MPOVA1.VPOCHA(/1)
      SEGACT MELEP1
 
      IKAS=1
      CALL XSOUR(FN,IZF1.FN,GR,PG,XYZ,HR,PGSQ,RPG,
     & NES,IDIM,NP,MP1,NPG,IAXI,IPT1.NUM,IKAS,KPRE,
     & IZTGG3.VPOCHA,IKS,NELG,IPADQ.LECT,MELEP1.NUM,
     & IZTGG3.VPOCHA,IKS,IZTGG3.VPOCHA,IKS,IPADS.LECT,
     & NBEL,NUTOEL,XCOOR,MPOVA1.VPOCHA,NPT)
 
      NUTOEL=NUTOEL+NBEL
      SEGSUP IZFFM,IZHR,IZF1,IZH2
      SEGDES IPT1
 1102 CONTINUE
      SEGDES IZTGG3,MPOVA1,MELEP1
      SEGSUP IPADQ
      IF(NBSOU1.NE.1)SEGDES MELEMK
 
      SEGACT MELEME
      NBSOU1=LISOUS(/1)
      IF(NBSOU1.EQ.0)NBSOU1=1
      LL=0
 
      IF(KIMPL.EQ.0.OR.KIMPL.EQ.2)THEN
      SEGACT MPOVA1*MOD
      DO 1533 L=1,NBSOU1
      IPT1=MELEME
      IF(NBSOU1.NE.1)IPT1=LISOUS(L)
      SEGACT IPT1
      NP=IPT1.NUM(/1)
      NBELEM=IPT1.NUM(/2)
 
      NNN=MOD(NBELEM,LRV)
      IF(NNN.EQ.0) NPACK=NBELEM/LRV
      IF(NNN.NE.0) NPACK=1+(NBELEM-NNN)/LRV
 
      DO 1534 KPACK=1,NPACK
 
C --- POUR CHAQUE PAQUET DE LRV ELEMENTS ou moins
      KDEB=1+(KPACK-1)*LRV
      KFIN=MIN(NBELEM,KDEB+LRV-1)
      NBEL=KFIN-KDEB+1
      LL=LL+1
 
      DO 2 N=1,NINKO
      IPMS=LIZAFS(LL,N)
      SEGACT IPMS
      DO 12 K=1,NBEL
      KI=KDEB+K-1
      DO 13 J=1,NP
      UU=0.D0
      IU=IPADS.LECT(IPT1.NUM(J,KI))
      DO 14 I=1,NP
      IK=IPADI.LECT(IPT1.NUM(I,KI))
      UU=UU+IPMS.AM(K,I,J)*IZTU1.VPOCHA(IK,N)
 14   CONTINUE
      MPOVA1.VPOCHA(IU,N)=MPOVA1.VPOCHA(IU,N)+UU
 13   CONTINUE
 12   CONTINUE
      SEGSUP IPMS
 2    CONTINUE
 
 1534 CONTINUE
      SEGDES IPT1
 1533 CONTINUE
      IF(NBSOU1.NE.1)SEGDES MELEME
      ELSE
      NBSOUS=LIZAFS(/1)
      DO      L=1,NBSOUS
      DO      N=1,NINKO
      IPMS=LIZAFS(L,N)
      SEGSUP IPMS
      ENDDO
      ENDDO
      ENDIF
      SEGDES IZTU1,MPOVA1
 
      CALL ECROBJ('MATRIK',MATRIK)
 
      CALL ECROBJ('CHPOINT',MCHPO1)
 
      SEGDES IMATRI
      SEGSUP IMATRS
      SEGDES MELEME,MATRIK
      IF(IKL.EQ.0)THEN
      SEGDES VISCO
      ENDIF
 
C*************************************************************************
      ELSE IF(KFORM.EQ.2)THEN
C CAS FORMULATION VF
      WRITE(6,*)' FORMULATION VF NON DISPONIBLE '
C        Option %m1:8 incompatible avec les données
            MOTERR( 1: 8) = '  EF    '
            CALL ERREUR(803)
            RETURN
 
      ENDIF
 
C     write(6,*)' RETOUR DE TSCA '
      IF(IPADS.NE.IPADI)THEN
      SEGSUP IPADS,IPADI
      ELSE
      SEGSUP IPADS
      ENDIF
      IPAS=1
      RETURN
 
 
 90   CONTINUE
      WRITE(6,*)' Interuption anormale de TSCAL '
C        Option %m1:8 incompatible avec les données
            CALL ERREUR(803)
      RETURN
 1002 FORMAT(10(1X,1PE11.4))
 1001 FORMAT(20(1X,I5))
      END
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

© Cast3M 2003 - Tous droits réservés.
Mentions légales