Télécharger zns.eso

Retour à la liste

zns
C ZNS       SOURCE    FANDEUR   22/01/03    21:16:04     11136          
      SUBROUTINE ZNS(MTABX,MTAB1)
      PARAMETER (NTB=1)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
C***********************************************************************
C
C  CET OPERATEUR DISCRETISE LES EQUATIONS DE NAVIER STOKES
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 du/dt + u Grad u = nu Lapl u - 1/ro Grad p  < + S >
C                                             < + g beta (T-Tref) >
C
C   'OPER' 'NS' nu            'INCO' UN :
C   'OPER' 'NS' nu s          'INCO' UN :
C   'OPER' 'NS' nu gb tn tref 'INCO' UN :
C
C
C    2/ Cas compréssible
C
C dG/dt + u Grad G + G Div u = mu Lapl u - Grad p   < + S >
C
C   'OPER' 'NS' mu un    'INCO' GN :
C   'OPER' 'NS' mu un s  'INCO' GN :
C
C
C    nu,mu viscosité cinématique resp. dynamique
C         FLOTTANT ou CHPOINT SCAL CENTRE
C    s    source volumique de qdm
C         POINT ou CHPOINT VECT CENTRE
C    gb   coéfficient de flottabilité (g*beta où g est l'accéllération
C         de la pesanteur et beta le coéfficient de dilatabilité)
C         POINT ou CHPOINT VECT CENTRE
C    tn   Champ de température CHPOINT SCAL SOMMET
C    tref température de référence
C         FLOTTANT ou CHPOINT SCAL SOMMET
C
C         Champ de vitesse                       -> VITESS
C    un   Champ de vitesse transportant          -> UTRANS
C         CHPOINT VECT SOMMET
C    gn   Champ de vitesse massique (transporté) -> IZTU1  (Inconnue)
C         CHPOINT VECT SOMMET
C
C************************************************************************
 
C?-INC CCVQUA4
 
-INC PPARAM
-INC CCOPTIO
-INC CCGEOME
-INC SIZFFB
      POINTEUR IZF1.IZFFM
-INC SMCHAML
-INC SMCOORD
-INC SMLENTI
      POINTEUR IPADI.MLENTI,IPADU.MLENTI,IPADF.MLENTI,IPADS.MLENTI
      POINTEUR IPADQ.MLENTI
-INC SMELEME
      POINTEUR MELEM1.MELEME,MELEMC.MELEME,IGEOM0.MELEME,MELEMS.MELEME
      POINTEUR MELEMI.MELEME,MELEP1.MELEME
-INC SMCHPOI
      POINTEUR IZTU1.MPOVAL,IZGG1.MPOVAL
      POINTEUR VITESS.MPOVAL,UTRANS.MPOVAL
      POINTEUR VISCO.MPOVAL,IZTGG2.MPOVAL
      POINTEUR IZTGG3.MPOVAL,IZTGG4.MPOVAL
      POINTEUR IZVOL.MPOVAL,IZTCO.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,TYPE,TYPC,NOM0,NOMA,NOM
      CHARACTER*8 LTAB(NTB)
      DIMENSION KTAB(NTB),IXV(3),RO(1)
      SAVE IPAS
      DATA LTAB/'KIZX    '/,IPAS/0/,RO/1.D0/
C*****************************************************************************
CNS
C     write(6,*)' DEBUT NS '
 
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*****************************************************************************
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
      CALL ACME(KOPTI,'IDCEN',IDCEN)
      CALL ACME(KOPTI,'IKOMP',IKOMP)
      CALL ACME(KOPTI,'ALE',KKALE)
      CALL ACME(KOPTI,'KIMPL',KIMPL)
      CALL ACME(KOPTI,'KPOIN',KPRE)
      CALL ACME(KOPTI,'KFORM',KFORM)
 
      IF(KFORM.NE.0.AND.KFORM.NE.1)THEN
C        Option %m1:8 incompatible avec les données
            MOTERR( 1: 8) = 'EF/EFM1 '
            CALL ERREUR(803)
            RETURN
        ENDIF
      CALL ACME(KOPTI,'KMACO',KMACO)
      CALL ACMF(KOPTI,'AIMPL',AIMPL)
      IF (IERR.NE.0) RETURN
 
C     write(6,*)' Apres les options '
C*****************************************************************************
C
C- Récupération de la table DOMAINE associée au domaine local
C
      CALL ACMM(MTABX,'NOMZONE',NOMZ)
      CALL LEKTAB(MTABX,'DOMZ',MTABZ)
      IF(MTABZ.EQ.0)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
 
      CALL LEKTAB(MTABZ,'MAILLAGE',MELEME)
      CALL LEKTAB(MTABZ,'SOMMET',MELEMS)
      CALL LEKTAB(MTABZ,'CENTRE',MELEMC)
      CALL LEKTAB(MTABZ,'MACRO',MACRO)
      CALL LEKTAB(MTABZ,'MACRO1',MELEMI)
      CALL LEKTAB(MTABZ,'QUADRATI',MQUAD)
      IF (IERR.NE.0) RETURN
 
      MELEMQ=MELEMC
      MELEP1=MELEMC
C     write(6,*)' KPRE=',kpre,' MACRO=',macro,' QUADR=',MQUAD
C?    IF(KPRE.LE.3.AND.MACRO.EQ.0.AND.MQUAD.EQ.0.AND.KFORM.NE.0)THEN
C?    WRITE(6,*)' Operateur NS '
C?    WRITE(6,*)' Incompatibilité du terme source et des éléments'
C?    WRITE(6,*)' MACRO ou QUADRATIQUE attendu '
C?    WRITE(6,*)' Interuption anormale de NS'
C        Option %m1:8 incompatible avec les données
C?          MOTERR( 1: 8) = '  EF    '
C?    CALL ERREUR(803)
C?    RETURN
          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
 
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
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 LICHT(MCHPOI,IZTU1,TYPC,MELEM1)
         NINKO = IZTU1.VPOCHA(/2)
         IF (NINKO.NE.IDIM) 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
C     write(6,*)' MCHPOI,MELEM1=',MCHPOI,MELEM1
C On fait pointer ces deux tableaux sur le champ U inconu (tjs présent) pour
C eviter de les enlever lors de l'appel FORTRAN si les options sont absentes
      UTRANS=IZTU1
      IKW=0
      VITESS=IZTU1
      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(MELEM1,IPADI)
      IPADU=IPADI
      IPADF=IPADI
      IPADS=IPADI
      IF(MELEM1.NE.MELEMS)CALL KRIPAD(MELEMS,IPADS)
 
      IF(IPAS.EQ.0)THEN
      IPAS=1
      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)
         CALL ERREUR(788)
      RETURN
      ENDIF
      ENDIF
 
C*****************************************************************************
 
 
 
C*************************************************************************
C Lecture des coefficient
C Type du coefficient :
C IK1=0 CHPOINT IK1=1 scalaire  IK1=2 vecteur
 
      CALL ACME(MTABX,'IARG',IARG)
      IF(IKOMP.EQ.0)THEN
      IF(IARG.NE.1.AND.IARG.NE.2.AND.IARG.NE.3.AND.IARG.NE.4.
     & AND.IARG.NE.5)THEN
      WRITE(6,*)' Operateur NS : option incompréssible '
      WRITE(6,*)' Nombre d''arguments incorrect : ',IARG
      WRITE(6,*)' On attend 1 , 2 ou 4 '
C           Indice %m1:8 : nombre d'arguments incorrect
            MOTERR(1:8) = 'IARG    '
            CALL ERREUR(804)
      RETURN
      ENDIF
      ELSEIF(IKOMP.EQ.1)THEN
      IF(IARG.NE.2.AND.IARG.NE.3)THEN
      WRITE(6,*)' Operateur NS : option compréssible '
      WRITE(6,*)' Nombre d''arguments incorrect : ',IARG
      WRITE(6,*)' On attend 2 ou 3 '
C           Indice %m1:8 : nombre d'arguments incorrect
            MOTERR(1:8) = 'IARG    '
            CALL ERREUR(804)
      RETURN
      ENDIF
      ENDIF
 
      IF(KKALE.NE.0)IARG=IARG-1
C--Cas incompréssible
      IF(IKOMP.EQ.0)THEN
      IXV(1)=MELEMC
      IXV(2)=1
      IXV(3)=0
      CALL LEKCOF('Opérateur NS :',
     & MTABX,KINC,1,IXV,IZTG1,VISCO,NPT1,NC1,IK1,IRET)
      IF(IRET.EQ.0)RETURN
 
      IK2=-1
      IKG=-1
      IK3=-1
      IK4=-1
      IZTGG2=VISCO
      IZTGG3=VISCO
      IZTGG4=VISCO
 
C 2ème coef -> source S
      IF(IARG.GE.2)THEN
      IXV(1)=-MELEMQ
      IXV(2)=0
      IXV(3)=1
      CALL LEKCOF('Opérateur NS :',
     & MTABX,KINC,2,IXV,IZTG2,IZTGG2,NELG,NC2,IKG,IRET)
      IF(IRET.EQ.0)RETURN
 
      IF(IKG.EQ.2)IKG=1
 
 
      IF(IARG.EQ.4)THEN
 
      IXV(1)=MELEMS
      IXV(2)=0
      IXV(3)=0
      CALL LEKCOF('Opérateur NS :',
     & MTABX,KINC,3,IXV,IZTG3,IZTGG3,NPT3,NC3,IK3,IRET)
      IF(IRET.EQ.0)RETURN
 
      IXV(1)=MELEMS
      IXV(2)=1
      IXV(3)=0
      CALL LEKCOF('Opérateur NS :',
     & MTABX,KINC,4,IXV,IZTG4,IZTGG4,NPT4,NC4,IK4,IRET)
      IF(IRET.EQ.0)RETURN
 
      ENDIF
      ENDIF
 
C--Cas compréssible
      ELSEIF(IKOMP.EQ.1)THEN
C 1er coef : mu , viscosité dynamique
      IXV(1)=MELEMC
      IXV(2)=1
      IXV(3)=0
      CALL LEKCOF('Opérateur NS :',
     & MTABX,KINC,1,IXV,IZTG1,VISCO,NPT1,NC1,IK1,IRET)
      IF(IRET.EQ.0)RETURN
 
      IK2=-1
      IKG=-1
      IK3=-1
      IK4=-1
      IZTGG2=VISCO
      IZTGG3=VISCO
      IZTGG4=VISCO
 
      IF(IARG.EQ.3)THEN
C 3ème coef : rog
      IXV(1)=-MELEMC
      IXV(2)=0
      IXV(3)=1
      CALL LEKCOF('Opérateur NS :',
     & MTABX,KINC,3,IXV,IZTG2,IZTGG2,NELG,NC2,IKG,IRET)
      IF(IRET.EQ.0)RETURN
 
      IF(IKG.EQ.2)IKG=1
 
      ENDIF
 
C 2ème coef : un , champ de vitesse transportant
      IXV(1)=-MELEMS
      IXV(2)=0
      IXV(3)=0
      CALL LEKCOF('Opérateur NS :',
     & MTABX,KINC,2,IXV,MUTRAN,UTRANS,NPTU,NC3,IKW,IRET)
      IF(IRET.EQ.0)RETURN
      VITESS=UTRANS
      IPADU=IPADS
      IPADF=IPADS
 
      ENDIF
 
      IF(KKALE.NE.0)THEN
 
      IF(IKOMP.EQ.1)THEN
      WRITE(6,*)' Operateur NS : option compréssible en ALE'
      WRITE(6,*)' Cas non prévu pour l''instant: '
C        Option %m1:8 incompatible avec les données
            MOTERR( 1: 8) = '  EF    '
      CALL ERREUR(803)
      RETURN
      ENDIF
 
C n+ unième coefficient champ de vitesse transportant pour l'option ALE
      IXV(1)=-MELEMS
      IXV(2)=0
      IXV(3)=0
      CALL LEKCOF('Opérateur NS :',
     & MTABX,KINC,(IARG+1),IXV,MUTRAN,UTRANS,NPTW,NCW,IKW,IRET)
      IF(IRET.EQ.0)RETURN
      IPADU=IPADS
 
      ENDIF
 
C     write(6,*)' Operateur NS : Fin lecture Arguments '
C Fin lecture Arguments ************************************************
 
      IF(IPADS.EQ.0)CALL KRIPAD(MELEMS,IPADS)
 
      CALL LEKTAB(MTABZ,'XXDXDY',MCHPOI)
      IF(MCHPOI.EQ.0)GO TO 90
      CALL LICHT(MCHPOI,IZTCO,TYPC,IGEOM)
      NELZ=IZTCO.VPOCHA(/1)
 
      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.E30
      DTP=1.E30+DT
      IPAT=1
      CALL ECME(MTABT,'NUPASDT',IPAT)
      ELSE
      CALL ACMF(MTABT,'DELTAT',DTP)
      ENDIF
 
C*********** FORMULATIONS **********
 
      IF(KFORM.EQ.0)THEN
C Formulation EFM1
      IF(KIMPL.NE.0)THEN
      WRITE(6,*)' Operateur NS '
C        Option %m1:8 incompatible avec les données
            MOTERR( 1: 8) = '  EF    '
      CALL ERREUR(803)
      RETURN
      ENDIF
 
      CALL LEKTAB(MTABZ,'MATESI',MATRIK)
      IF(MATRIK.EQ.0)GO TO 90
      SEGACT MATRIK
 
      IMATRI=IRIGEL(4,1)
      SEGACT IMATRI
 
      CALL LEKTAB(MTABZ,'XXPSOML',MCHELM)
      IF(MCHELM.EQ.0)GO TO 90
      SEGACT MCHELM
 
      CALL LEKTAB(MTABZ,'XXVOLUM',MCHVOL)
      IF(MCHVOL.EQ.0)GO TO 90
      CALL LICHT(MCHVOL,IZVOL,TYPC,IGEOM)
 
      CALL LEKTAB(MTAB1,'KIZG',KIZG)
      IF(KIZG.EQ.0)THEN
      CALL CRTABL(KIZG)
      CALL ECMM(KIZG,'SOUSTYPE','KIZG')
      CALL ECMO(MTAB1,'KIZG','TABLE   ',KIZG)
      ENDIF
 
      TYPE=' '
      CALL ACMO(KIZG,NOMI,TYPE,IZG1)
      IF(TYPE.NE.'CHPOINT ')THEN
      NC=IZTU1.VPOCHA(/2)
      TYPE='SOMMET'
      CALL CRCHPT(TYPE,MELEM1,NC,IZG1)
      CALL ECMO(KIZG,NOMI,'CHPOINT ',IZG1)
      ENDIF
 
      CALL LICHT(IZG1,IZGG1,TYPC,IGEOM)
 
      IF(IGEOM.NE.MELEM1)THEN
      WRITE(6,*)' Opérateur NS'
      WRITE(6,*)'Incompatibilite de SPG entre 1eres inconnues'
C        Option %m1:8 incompatible avec les données
            MOTERR( 1: 8) = '  EF    '
      CALL ERREUR(803)
      RETURN
      ENDIF
 
      NPT=IZGG1.VPOCHA(/1)
 
      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
      IZAFM=LIZAFM(L,1)
      IPM1=IZAFM
      SEGACT IZAFM
      IF(IAXI.NE.0)THEN
      IPM1=LIZAFM(L,2)
      SEGACT IPM1
      ENDIF
 
      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
 
      NP  =IPT1.NUM(/1)
      NBEL=IPT1.NUM(/2)
      IES=IDIM
 
       NPTF=VITESS.VPOCHA(/1)
       NPTU=UTRANS.VPOCHA(/1)
       CALL ZCVI(AM,IPM1.AM,VELCHE,
     & IPT1.NUM,NBEL,NUTOEL,NPT,IES,NP,IAXI,IKOMP,IARG,
     & VISCO.VPOCHA,IK1,IZTGG2.VPOCHA,IKG,NELG,IZTGG3.VPOCHA,IK3,
     & IZTGG4.VPOCHA,IK4,IPADS.LECT,
     & UTRANS.VPOCHA,IPADU.LECT,NPTU,
     & IZTU1.VPOCHA,IZGG1.VPOCHA,IPADI.LECT,
     & VITESS.VPOCHA,IPADF.LECT,NPTF,
     & IZVOL.VPOCHA,IZTCO.VPOCHA,NELZ,IDCEN,IPG,
     & DTM1,DT,DTT1,DTT2,NUEL,DIAEL,IZFFM.FN)
 
      SEGDES IZAFM,IPT1
      IF(IAXI.NE.0)SEGDES IPM1
      NUTOEL=NUTOEL+NBEL
 
 1    CONTINUE
      SEGDES MELEME
 
 
      IF(DT.LT.DTP)THEN
      CALL ECMF(MTABT,'DELTAT',DT)
      CALL ECMM(MTABT,'OPER','NS')
      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,IZTGG2,IZTGG3,IZTGG4
      SEGDES UTRANS,IZTU1
      SEGDES IZGG1
      SEGDES IZVOL,IZTCO
      SEGDES LINCO
      SEGDES MATRIK,IMATRI
 
 
C*************************************************************************
      ELSE IF(KFORM.EQ.1)THEN
C CAS FORMULATION EF
 
      IF(KIMPL.EQ.0)THEN
C        Option %m1:8 incompatible avec les données
            MOTERR( 1: 8) = '  EF    '
            CALL ERREUR(803)
            RETURN
        ENDIF
 
      DT=0.D0
      IF(IDCEN.EQ.4)THEN
      TYPE=' '
      CALL ACMO(MTABT,'DELTAT',TYPE,IENT)
      IF(TYPE.NE.'ENTIER')THEN
      WRITE(6,*)' Opérateur NS '
      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=1
      NUTOEL=0
      NINKO=VITESS.VPOCHA(/2)
      SEGACT MELEME
      NBSOUS=LISOUS(/1)
      IF(NBSOUS.EQ.0)NBSOUS=1
 
      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
      NBOP=0
      NBME=NINKO
      NBELC=0
      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
 
      NUTOEL=0
      DO 101 L=1,NBSOUS
      IPT1=MELEME
      IF(NBSOUS.NE.1)IPT1=LISOUS(L)
      SEGACT IPT1
      NOM0=NOMS(IPT1.ITYPEL)//'    '
      CALL KALPBG(NOM0,'FONFORM ',IZFFM)
      SEGACT IZFFM*MOD
      IZHR=KZHR(1)
      SEGACT IZHR*MOD
      NES=GR(/1)
      NPG=GR(/3)
 
      NP = IPT1.NUM(/1)
      MP = NP
      NBEL=IPT1.NUM(/2)
      SEGINI IPM1,IPS1
      LIZAFM(L,1)=IPM1
      LIZAFS(L,1)=IPS1
      IPM2=IPM1
      IPM3=IPM1
      IPS2=IPS1
      IPS3=IPS1
      IF(NBME.GE.2)THEN
      SEGINI IPM2,IPS2
      LIZAFM(L,2)=IPM2
      LIZAFS(L,2)=IPS2
      ENDIF
      IF(NBME.GE.3)THEN
      SEGINI IPM3,IPS3
      LIZAFM(L,3)=IPM3
      LIZAFS(L,3)=IPS3
      ENDIF
 
 
      KITT=2
      KJTT=IK1
      NPT=UTRANS.VPOCHA(/1)
      NPTG=IZTU1.VPOCHA(/1)
 
      CALL ZCONV(FN,GR,PG,XYZ,HR,PGSQ,RPG,
     & NES,IDIM,NP,NPG,IAXI,AIMPL,IKOMP,
     & RO,1,UTRANS.VPOCHA,IKW,NPT,IPADU.LECT,VISCO.VPOCHA,IK1,
     & IPT1.NUM,NBEL,NUTOEL,XCOOR,
     & IPM1.AM,IPM2.AM,IPM3.AM,
     & IPS1.AM,IPS2.AM,IPS3.AM,
     & NINKO,IDCEN,DT,
     & IZTU1.VPOCHA,0,NPTG,IPADI.LECT,IZTCO.VPOCHA,NELZ)
 
      CALL XLAPL(FN,GR,PG,XYZ,HR,PGSQ,RPG,NES,IDIM,NP,NPG,IAXI,
     & VISCO.VPOCHA,VISCO.VPOCHA,VISCO.VPOCHA,KITT,KJTT,IK1,
     & IPT1.NUM,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
 101  CONTINUE
 
      ENDIF
 
 
      IF(KIMPL.EQ.2.OR.KIMPL.EQ.0.OR.IARG.GT.1)THEN
      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
      ENDIF
 
      IF(IARG.EQ.2.OR.IARG.EQ.4)THEN
       CALL KRIPAD(MELEMQ,IPADQ)
      IF(IARG.EQ.2)THEN
       IKAS=2
      ELSEIF(IARG.EQ.4)THEN
       IKAS=3
      ENDIF
 
 
 
      IF(MACRO.NE.0.AND.KPRE.EQ.2)MELEMI=MELEME
      SEGACT MELEMI
      NBSOUS=MELEMI.LISOUS(/1)
      IF(NBSOUS.EQ.0)NBSOUS=1
 
      NUTOEL=0
      DO 1102 L=1,NBSOUS
      IPT1=MELEMI
      IF(NBSOUS.NE.1)IPT1=MELEMI.LISOUS(L)
      SEGACT IPT1
 
       IF(MQUAD.NE.0)THEN
       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(MACRO.NE.0)THEN
       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)
      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=IZTGG2.VPOCHA(/1)
      NPT =MPOVA1.VPOCHA(/1)
 
      SEGACT MELEP1
 
 
      CALL XSOUR(FN,IZF1.FN,GR,PG,XYZ,HR,PGSQ,RPG,
     & NES,IDIM,NP,MP1,NPG,IAXI,IPT1.NUM,IKAS,KPRE,
     & IZTGG2.VPOCHA,IKG,NELG,IPADQ.LECT,MELEP1.NUM,
     & IZTGG3.VPOCHA,IK3,IZTGG4.VPOCHA,IK4,IPADS.LECT,
     & NBEL,NUTOEL,XCOOR,MPOVA1.VPOCHA,NPT)
 
      NUTOEL=NUTOEL+NBEL
 1102 CONTINUE
 
      ENDIF
 
      IF(KIMPL.EQ.2.OR.KIMPL.EQ.0)THEN
 
      NBSOUS=LISOUS(/1)
      IF(NBSOUS.EQ.0)NBSOUS=1
 
      DO 1533 L=1,NBSOUS
      IPT1=MELEME
      IF(NBSOUS.NE.1)IPT1=LISOUS(L)
      SEGACT IPT1
      NP=IPT1.NUM(/1)
      NBEL=IPT1.NUM(/2)
      DO 2 N=1,NINKO
      IPMS=LIZAFS(L,N)
      SEGACT IPMS
      DO 12 K=1,NBEL
      DO 13 J=1,NP
      UU=0.D0
      IU=IPADS.LECT(IPT1.NUM(J,K))
      DO 14 I=1,NP
      IK=IPADI.LECT(IPT1.NUM(I,K))
      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
 
 2    CONTINUE
 
 1533 CONTINUE
      ENDIF
 
      SEGDES IPM1,IPM2,IPM3
      IPS=IPS1
      SEGSUP IPS1
      IF(IPS2.NE.IPS)SEGSUP IPS2
      IF(IPS3.NE.IPS)SEGSUP IPS3
      SEGDES IZTCO
 
      IF(KIMPL.EQ.2.OR.KIMPL.EQ.0.OR.IARG.GT.1)THEN
 
      TYPE=' '
      CALL ACMO(MTAB1,'SMBR',TYPE,MCHPO2)
      IF(TYPE.NE.'CHPOINT')THEN
      CALL ECMO(MTAB1,'SMBR','CHPOINT',MCHPO1)
      ELSE
      CALL ECROBJ('CHPOINT',MCHPO2)
      CALL ECROBJ('CHPOINT',MCHPO1)
C?    CALL OPERAD
      CALL PRFUSE
      CALL LIROBJ('CHPOINT',MCHPOI,1,IRET)
      CALL ECMO(MTAB1,'SMBR','CHPOINT',MCHPOI)
      ENDIF
 
      ENDIF
 
      SEGDES IMATRI
      SEGDES MELEME,MATRIK
      IF(IK1.EQ.0)THEN
      SEGDES VISCO
      ENDIF
      CALL ECMO(MTABX,'MATELM','MATRIK',MATRIK)
 
C*************************************************************************
      ELSE IF(KFORM.EQ.2)THEN
C CAS FORMULATION VF
      WRITE(6,*)' FORMULATION VF NON DISPONIBLE '
      ENDIF
C*************************************************************************
 
      IPDI=IPADI
      SEGSUP IPADI
      IF(IPADS.NE.IPDI)SEGSUP IPADS
      IPAS=1
C     write(6,*)' FIN de NS'
      RETURN
 90   CONTINUE
      WRITE(6,*)' Interuption anormale de NS'
C        Option %m1:8 incompatible avec les données
            MOTERR( 1: 8) = '  EF    '
      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