Télécharger modeli.eso

Retour à la liste

modeli
C MODELI    SOURCE    OF166741  24/12/18    21:15:25     12090          
C----------------------------------------------------------------------C
C                     OPERATEUR MODELE                                 C
C                                                                      C
C  Creation d'un objet MODELE                                          C
C                                                                      C
C  Syntaxe : MOD1 = MODL GEO1 TYPE_CAL TYPE_MAT ( TYPE_ELE ) ;         C
C                                                                      C
C    GEO1        MAILLAGE de base                                      C
C    TYPE_CAL    MOT(S) pour definir la FORMULATION                    C
C    TYPE_MAT    MOT(S) pour definir le MATERIAU                       C
C    TYPE_ELE    MOT(S) pour definir les ELEMENTS FINIS a utiliser     C
C    MOD1        Resultat de type MODELE                               C
C----------------------------------------------------------------------C
C  PPU : Modif pour les materiaux unidirectionels en plasticite
 
      SUBROUTINE MODELI
 
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
 
-INC PPARAM
-INC CCOPTIO
-INC CCHAMP
-INC CCGEOME
C==DEB= FORMULATION HHO == Include specifique ==========================
-INC CCHHOPA
C==FIN= FORMULATION HHO ================================================
 
-INC SMCOORD
-INC SMELEME
-INC SMMODEL
      POINTEUR IMODE3.IMODEL,IMODE4.IMODEL,IMODE5.IMODEL
      POINTEUR nomid1.NOMID,nomid2.NOMID
-INC SMTABLE
-INC SMLMOTS
      POINTEUR OPNLIN.MLMOTS
 
      SEGMENT LIMODE(0)
      SEGMENT PLICON
        integer mlicon(NLCON),tlicon(NLCON)
      ENDSEGMENT
 
      EXTERNAL LONG
      PARAMETER (NBFORM=19,NBCON=13,NBEXT=7,NBDIF=1)
      PARAMETER (N1MAX=300,N2MAX=200)
      PARAMETER (NLOMAX=5)
 
      DIMENSION LESMOD(N1MAX)
      CHARACTER*4 MOTEF(N2MAX),LESTEF(N2MAX),MOCON(NBCON),MOEXT(NBEXT),
     &            MOINCO(NBDIF)
      CHARACTER*4 MNLOCA(NLOMAX),MNLVAR(1)
      CHARACTER*4 MCTCT(4)
      CHARACTER*8 TAPIND,TYPOBJ,CHARIN,CHARRE,CMATE,PHAM
      CHARACTER*8 PAR1,MDIINC,MDIDUA
      CHARACTER*(LCONMO) CONM
      CHARACTER*(LOCOMP) MOPRID
      CHARACTER*16 MOFORM(NBFORM),LESFOR(2),MOPROP(N1MAX),LESPRO(N1MAX)
      CHARACTER*16 LMENOM,LDINOM,OPTEMP(3)
      CHARACTER*(LOCHAI) MOTEMP,LMELIB,LDILIB,LMEFCT,LDIFCT
      CHARACTER*4 mgauss(4)
      CHARACTER*4 deriv(1)
 
      LOGICAL LOGRE,LOGIN,LMEEXT,LMENLX,LMEVIX,LOSTAT,LOMELA,LINOMID
      LOGICAL LDIEXT,LDISOR,LOBBAR
 
C=DEB==== FORMULATION HHO ==== Declarations particulieres ==============
      PARAMETER (NMHHO=2)
      CHARACTER*4 mcHHO(NMHHO)
      CHARACTER*(LOCHAI) chaHHO
      LOGICAL loHHO
      DATA mcHHO / 'HHO_','HHO ' /
C=FIN==== FORMULATION HHO ==============================================
 
      DATA MGAUSS /'EPAI' , 'RIGI' , 'MASS' ,'CONT'/
      DATA DERIV  /'EPSI'/
c       DATA MDERIV/'LINEAIRE        ','QUADRATIQUE     ',
c      $            'TRUESDELL       ','JAUMANN         ',
c      $            'UTILISATEUR     ','FEFP            '/
      DATA OPTEMP/'PHASE           ','ADVECTION       ',
     $            'CONDUCTION      '/
 
C----------------------------------------------------------------------C
C              DEFINITION DES NOMS DE FORMULATIONS                     C
C  Formulation LIAISON : pour  operateurs DYNE et COMP                 C
C----------------------------------------------------------------------C
      DATA MOFORM /
     &     'THERMIQUE       ','MECANIQUE       ','LIQUIDE         ',
     &     'POREUX          ','DARCY           ','CONTACT         ',
     &     'MAGNETODYNAMIQUE','NAVIER_STOKES   ','MELANGE         ',
     &     'EULER           ','FISSURE         ','LIAISON         ',
     &     'THERMOHYDRIQUE  ','ELECTROSTATIQUE ','DIFFUSION       ',
     &     'CHARGEMENT      ','METALLURGIE     ','CHANGEMENT_PHASE',
     &     'CONTRAINTE      '                                      /
 
C (fdp) Ajout d'un nouveau mot clef 'LIBRE' ou 'LIE' pour les JOI1
      DATA MOCON / 'CONS','INTE','DPGE','PHAS','STAT','LCOI','LCOS',
     &             'LIBR','LIE ','NON_','LINE','CHPO','GAP7'/
      DATA MOEXT / 'NUME','NOM_','PARA','C_MA','C_VA','LIB_','FCT_' /
      DATA MOINCO / 'INCO' /
      DATA MNLVAR/ 'V_MO' /
      DATA MCTCT/'MESC','FAIB','SYME','MORT'/
 
      CONM  ='                '
      PHAM  ='        '
      MDIINC='        '
      MDIDUA='        '
 
      MFR    = 0
      lucvar = 0
      lucmat = 0
      lucmaf = 0
      luparx = 0
      lobbar = .false.
      lecont = 0
C=DEB==== FORMULATION HHO ==== Initialisations particulieres ===========
      loHHO  = .FALSE.
C=FIN==== FORMULATION HHO ==============================================
 
      mmode2 = 0
 
      IPGEOM = 0
      IPTABL = 0
      IPTABM = 0
 
C  Lecture d'une table STATIONNAIRE
      IPTABL = 0
      IPTABS = 0
      IPTABM = 0
      IPGEOM = 0
      CALL LIRTAB('STATIONNAIRE',IPTABL,0,IRET)
      IF (IERR.NE.0) RETURN
      IF (IRET.GT.0) THEN
        IPTABS=IPTABL
        IVALIN=0
        XVALIN=REAL(0.D0)
        LOGIN=.TRUE.
        IOBIN=0
        TAPIND='MOT     '
        CHARIN='MAILLAGE'
        TYPOBJ='TABLE   '
        CALL ACCTAB(IPTABS,TAPIND,IVALIN,XVALIN,CHARIN,LOGIN,IOBIN,
     .                     TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE)
        IF (IERR.NE.0) RETURN
        IPTABM = IOBRE
        IVALIN=1
        XVALIN=REAL(0.D0)
        LOGIN=.TRUE.
        IOBIN=0
        TAPIND='ENTIER  '
        CHARIN='        '
        TYPOBJ='MAILLAGE'
        CALL ACCTAB(IPTABM,TAPIND,IVALIN,XVALIN,CHARIN,LOGIN,IOBIN,
     .                     TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE)
        IF (IERR.NE.0) RETURN
        IPGEOM = IOBRE
        IRET = 0
      ENDIF
 
C Lecture d'une TABLE de sous-type MAILLAGE
      IF(IPTABM.EQ.0) THEN
       CALL LIRTAB('MAILLAGE',IPTABL,0,IRET)
       IF (IERR.NE.0) RETURN
       IF (IRET.GT.0) THEN
        IPTABM = IPTABL
        IVALIN=1
        XVALIN=REAL(0.D0)
        LOGIN=.TRUE.
        IOBIN=0
        TAPIND='ENTIER  '
        CHARIN='        '
        TYPOBJ='MAILLAGE'
        CALL ACCTAB(IPTABM,TAPIND,IVALIN,XVALIN,CHARIN,LOGIN,IOBIN,
     .                     TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE)
        IF (IERR.NE.0) RETURN
        IPGEOM = IOBRE
        IRET = 0
       ENDIF
      ENDIF
 
C  Lecture d'un MAILLAGE ou d'une TABLE de sous-type DOMAINE
      IPTABL = 0
      IPGEO2 = 0
      IReMOD = 0
      CALL LIRTAB('DOMAINE',IPTABL,0,IRET)
      IF (IERR.NE.0) RETURN
      IF (IPTABL.GT.0) THEN
        IVALIN=0
        XVALIN=REAL(0.D0)
        LOGIN=.TRUE.
        IOBIN=0
        TAPIND='MOT     '
        CHARIN='MAILLAGE'
        TYPOBJ='MAILLAGE'
        CALL ACCTAB(IPTABL,TAPIND,IVALIN,XVALIN,CHARIN,LOGIN,IOBIN,
     &                     TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE)
        IF (IERR.NE.0) RETURN
        IPGEOM = IOBRE
      ENDIF
 
C  Lecture d'un MAILLAGE (cas general) :
      IF (IPGEOM.LE.0) THEN
        CALL LIROBJ('MAILLAGE',IPGEOM,1,IRET)
        IF (IERR.NE.0) RETURN
      ENDIF
 
C  Verification de l'unicite des elements
c  On ne tient pas compte de l'ordre des noeuds dans l'element
      IPT1 = IPGEOM
      iordre=0
      CALL UNIQMA(IPT1,NBDI1,iordre)
      IF (NBDI1 .NE. 0) THEN
        MOTERR(1:8)='MAILLAGE'
        CALL ERREUR(1019)
        RETURN
      ENDIF
 
C  Lecture d'une FORMULATION
      NFOR =0
      NMAT =0
      CALL MESLIR(-182)
      ICOND=1
      CALL MESLIR(-182)
 
 51   CONTINUE
      CALL LIRMOT(MOFORM,NBFORM,IPFORM,ICOND)
      IF (IERR.NE.0) RETURN
      IF (IPFORM .EQ. 0) GOTO 52
 
      NFOR=NFOR+1
      IF (NFOR.GT.2) THEN
        CALL ERREUR(251)
        RETURN
      ENDIF
      LESFOR(NFOR)=MOFORM(IPFORM)
      ICOND=0
      CALL MESLIR(-181)
      GOTO 51
 
 52   CONTINUE
C  Cas d'une FORMULATION simple (NFOR=1)
      IF (NFOR.EQ.1) THEN
c         jderiv=mepsil
cbp,2020-12-10 : abandon de MEPSIL (CCOPTIO) et IDERIV (MMODEL)
        IF     (LESFOR(1).EQ.'THERMIQUE') THEN
          CALL MODEL1(MOPROP,NPROP,MOTEF,NBTEF,N1MAX,N2MAX)
        ELSE IF(LESFOR(1).EQ.'MECANIQUE') THEN
          CALL MODEL2(MOPROP,NPROP,MOTEF,NBTEF,N1MAX,N2MAX)
        ELSE IF(LESFOR(1).EQ.'LIQUIDE') THEN
          CALL MODEL3(MOPROP,NPROP,MOTEF,NBTEF,N1MAX,N2MAX)
        ELSE IF(LESFOR(1).EQ.'POREUX') THEN
          CALL MODEL6(MOPROP,NPROP,MOTEF,NBTEF,N1MAX,N2MAX)
        ELSE IF(LESFOR(1).EQ.'DARCY') THEN
          CALL MODEL7(MOPROP,NPROP,MOTEF,NBTEF,N1MAX,N2MAX)
        ELSE IF(LESFOR(1).EQ.'CONTACT') THEN
          CALL MODEL8(MOPROP,NPROP,MOTEF,NBTEF,N1MAX,N2MAX)
        ELSE IF(LESFOR(1).EQ.'MAGNETODYNAMIQUE') THEN
          CALL MODE10(MOPROP,NPROP,MOTEF,NBTEF,N1MAX,N2MAX)
        ELSE IF(LESFOR(1).EQ.'NAVIER_STOKES') THEN
          CALL MODE11(MOPROP,NPROP,MOTEF,NBTEF,N1MAX,N2MAX)
        ELSE IF (LESFOR(1).EQ.'MELANGE') THEN
          CALL MODE12(MOPROP,NPROP,MOTEF,NBTEF,N1MAX,N2MAX)
          DO i=1,N1MAX
            LESMOD(i)=0
          ENDDO
        ELSE IF(LESFOR(1).EQ.'EULER') THEN
          CALL MODE13(MOPROP,NPROP,NBTEF,N1MAX)
        ELSE IF(LESFOR(1).EQ.'FISSURE') THEN
          CALL MODE14(MOPROP,NPROP,MOTEF,NBTEF,N1MAX,N2MAX)
        ELSE IF(LESFOR(1).EQ.'LIAISON') THEN
          CALL MODE15(MOPROP,NPROP,MOTEF,NBTEF,N1MAX,N2MAX)
        ELSE IF(LESFOR(1).EQ.'THERMOHYDRIQUE') THEN
          CALL MODE16(MOPROP,NPROP,MOTEF,NBTEF,N1MAX,N2MAX)
        ELSE IF(LESFOR(1).EQ.'ELECTROSTATIQUE ') THEN
          CALL MODE17(MOPROP,NPROP,MOTEF,NBTEF,N1MAX,N2MAX)
        ELSE IF(LESFOR(1).EQ.'DIFFUSION       ') THEN
          CALL MODE18(MOPROP,NPROP,MOTEF,NBTEF,N1MAX,N2MAX)
        ELSE IF(LESFOR(1).EQ.'CHARGEMENT      ') THEN
          CALL MODE19(MOPROP,NPROP,MOTEF,NBTEF,N1MAX,N2MAX)
        ELSE IF(LESFOR(1).EQ.'METALLURGIE     ') THEN
cjk148537 : ce n'est pas l exemple a suivre
          CALL MODE21(MOPROP,NPROP,MOTEF,NBTEF,N1MAX,N2MAX)
        ELSE IF(LESFOR(1).EQ.'CHANGEMENT_PHASE') THEN
          CALL MODE22(MOPROP,NPROP,MOTEF,NBTEF,N1MAX,N2MAX)
        ELSE IF(LESFOR(1).EQ.'CONTRAINTE') THEN
          CALL MODE24(MOPROP,NPROP,MOTEF,NBTEF,N1MAX,N2MAX)
        ELSE
          CALL ERREUR (251)
        ENDIF
        IF (IERR.NE.0) RETURN
 
      ELSE
C  Cas d'une FORMULATION couplee (NFOR=2)
c         jderiv=mepsil
cbp,2020-12-10 : abandon de MEPSIL (CCOPTIO) et IDERIV (MMODEL)
        IF ((LESFOR(1).EQ.'LIQUIDE'.AND.LESFOR(2).EQ.'MECANIQUE').OR.
     &      (LESFOR(2).EQ.'LIQUIDE'.AND.LESFOR(1).EQ.'MECANIQUE')) THEN
          CALL MODEL5(NPROP,MOTEF,NBTEF,N2MAX)
        ELSE
          CALL ERREUR(251)
        ENDIF
        IF (IERR.NE.0) RETURN
      ENDIF
 
C  Lecture eventuelle des proprietes du MODELE de MATERIAU
      ifrtt  = 0
      IFROCA = 0
      ifacaf = 0
      isyme  = 0
      nbga   = 10
      nbdang = 3
      icavit = 0
      kjh    = 0
      ISRCE  = 0
      IRAYE  = 0
      ICONV  = 0
      NMAT   = 0
      iprop  = 0
      ipgeo2 = 0
      if (iptabm.gt.0.and.iptabs.eq.0) goto 674
      IF (NPROP.EQ.0) GOTO 43
      CALL MESLIR(-180)
 41   CONTINUE
      IF (NMAT .NE.0) CALL MESLIR(-179)
      CALL LIRMOT(MOPROP,NPROP,LAPROP,0)
      IF (IERR.NE.0) RETURN
 
C   ---------- Cas d'un MODELE de CONTACT
      IF(LESFOR(1).EQ.'CONTACT') then
*  si FROTTANT lecont=3
        if(laprop.eq.3) lecont=laprop
        if(laprop.eq.5) then
          ifrtt=0
          ifroca=1
C         call lirobj('MMODEL',IFROCA,1,iOK)
C         IF(ierr.NE.0) return
          Call lirobj('MAILLAGE',IBETON,1,IOK)
          IF(ierr.NE.0) return
        endif
        if(laprop.eq.4) then
          ifrtt=1
        endif
      ENDIF
 
C   ---------- Cas d'un MODELE de CONTRAINTE
      IF(LESFOR(1).EQ.'CONTRAINTE') then
        if(laprop.ne.0) then
          lactr=laprop
          NMAT=NMAT+1
          LESPRO(NMAT)=MOPROP(LAPROP)
        endif
        if (lactr.ne.0) goto 42
      ENDIF
 
C   ---------- Cas d'un MODELE de METALLURGIE
C              modele cree par T.L. en mai 2018
      IF (lesfor(1).eq.'METALLURGIE' .AND. NMAT.le.4 ) THEN
        NMAT=NMAT+1
        IF( laprop .eq. 1 ) THEN
C        On vas lire le LISTMOTS qui suit le mot-clef MOPROP(1)='PHASES'
          CALL LIROBJ('LISTMOTS', lucvar, 1, IRETOU)
          MLMOTS = lucvar
          segact MLMOTS
          NB_PHA = MLMOTS.MOTS(/2)
C         On remplira ensuite MATMOD() avec lespro()
cjk148537          lespro(laprop) = MOPROP(laprop)
 
        ELSEIF( laprop .eq. 2 ) THEN
C        On vas lire le LISTMOTS qui suit le mot-clef MOPROP(2)='REACTIFS'
          CALL LIROBJ('LISTMOTS', ireact, 1, IRETOU)
          MLMOT1 = ireact
          segact MLMOT1
          NB_REA = MLMOT1.MOTS(/2)
C         On remplira ensuite MATMOD() avec lespro()
cjk148537          lespro(laprop) = MOPROP(laprop)
 
        ELSEIF( laprop .eq. 3 ) THEN
C        On vas lire le LISTMOTS qui suit le mot-clef MOPROP(3)='PRODUITS'
          CALL LIROBJ('LISTMOTS', iprodu, 1, IRETOU)
          MLMOT2 = iprodu
          segact MLMOT2
          NB_PRO = MLMOT2.MOTS(/2)
C         On remplira ensuite MATMOD() avec lespro()
cjk148537          lespro(laprop) = MOPROP(laprop)
 
        ELSEIF( laprop .eq. 4 ) THEN
C        On vas lire le LISTMOTS qui suit le mot-clef MOPROP(4)='TYPE'
          CALL LIROBJ('LISTMOTS', lucmat, 1, IRETOU)
          MLMOT3 = lucmat
          segact MLMOT3
          NB_TYP = MLMOT3.MOTS(/2)
C         On remplira ensuite MATMOD() avec lespro()
          do jj = 1,nb_typ
           lespro(jj) = mlmot3.mots(jj)
          enddo
 
        ELSE
          CALL ERREUR(5)
          RETURN
        ENDIF
 
C       Les pointeurs lucvar et lucmat sont ensuite transmis a
C       inomid pour remplir les NOMID de l'objet MMODEL
C       Les pointeurs lucvar, ireact, iprodu, lucmat seront mis
C       dans le tableau IVAMOD de l'objet IMODEL
 
        IF(NMAT .lt. 4) THEN
C         On n'a pas encore recuperer toutes les donnees
          go to 41
 
        ELSE IF(NMAT .eq. 4) THEN
 
C         On emet une erreur si les MLMOTS 'REACTIFS', 'PRODUITS' et
C         'TYPES' n'ont pas ete luts
          if(ireact .le. 0 .OR. iprodu .le. 0 .OR. lucmat .le. 0) then
            CALL ERREUR(21)
            RETURN
          endif
 
C         Autant de produits que de reactifs
          if( NB_PRO .ne. NB_REA ) then
            CALL ERREUR(1078)
            RETURN
          endif
 
C         On initialise le MLMOTS des PHASES si celui ci n'a pas ete lu
          icompt = 0
          if( lucvar .le. 0) then
            icompt = 1
            NB_PHA = NB_REA + NB_PRO
            JGN = LOCOMP
            JGM = NB_PHA
            SEGINI, MLMOTS
            lucvar = MLMOTS
C           On remplira ensuite MATMOD() avec lespro()
            lespro(1) = MOPROP(1)
          endif
 
C         On a recuperer toutes les donnees, on effectue quelques tests
          do ipha = 1, NB_PRO
 
C           Produits differents du reactif pour chaque reaction
            if( MLMOT1.MOTS(ipha) .eq. MLMOT2.MOTS(ipha) ) then
              MOTERR(1:4)=MLMOT1.MOTS(ipha)
              MOTERR(5:8)=MLMOT2.MOTS(ipha)
              CALL ERREUR(1075)
              RETURN
            endif
 
            irphas = 0
            ipphas = 0
            CALL PLACE(MLMOTS.MOTS, NB_PHA, irphas, MLMOT1.MOTS(ipha))
            CALL PLACE(MLMOTS.MOTS, NB_PHA, ipphas, MLMOT2.MOTS(ipha))
C             Si le nom du produit ou du reactif n'a pas ete lu dans le
C             MLMOTS des PHASES :
C               On le rajoute        si lucvar n'avait pas ete lu
C               On emet une erreur   sinon
            if(irphas .eq. 0) then
              if( icompt .ge. 1 ) then
                MLMOTS.MOTS(icompt) = MLMOT1.MOTS(ipha)
                icompt = icompt + 1
              else
                MOTERR(1:4)=MLMOT1.MOTS(ipha)
                CALL ERREUR(1080)
                RETURN
              endif
            endif
            if(ipphas .eq. 0) then
              if( icompt .ge. 1 ) then
                MLMOTS.MOTS(icompt) = MLMOT2.MOTS(ipha)
                icompt = icompt + 1
              else
                MOTERR(1:4)=MLMOT2.MOTS(ipha)
                CALL ERREUR(1080)
                RETURN
              endif
            endif
 
          enddo
C         On corrige la taille de MLMOTS :
          if( icompt .ge. 1 ) then
            JGM = icompt - 1
            JGN = MLMOTS.MOTS(/1)
            SEGADJ, MLMOTS
          endif
 
C         Un type de reaction definit pour chaque reaction
          if( NB_TYP .ne. NB_PRO ) then
            CALL ERREUR(1077)
            RETURN
          endif
 
          LAPROP = 0
 
          segact,MLMOTS*NOMOD, MLMOT1*NOMOD, MLMOT2*NOMOD, MLMOT3*NOMOD
        ENDIF
      ENDIF
 
C   ---------- Cas d'un MODELE de THERMIQUE CONVECTION ou RAYONNEMENT
      IF (lesfor(1).eq.'THERMIQUE' .AND. kjh.eq.0) then
        IF (moprop(laprop).eq.'CONVECTION') then
          ICONV=1
          nmat=nmat+1
          kjh=1
          lespro(nmat)=moprop(laprop)
          call model4(MOPROP,NPROP,MOTEF,NBTEF,N1MAX,N2MAX)
          go to 41
        ELSE IF (moprop(laprop).eq.'RAYONNEMENT') then
          IRAYE=1
          kjh=1
          nmat=nmat+1
          lespro(nmat)=moprop(laprop)
          call model9(MOPROP,NPROP,MOTEF,NBTEF,N1MAX,N2MAX)
          go to 41
        ELSE IF (moprop(laprop).eq.'SOURCE') then
          ISRCE=1
          kjh=1
          nmat=nmat+1
          lespro(nmat)=moprop(laprop)
          call mode23(MOPROP,NPROP,MOTEF,NBTEF,N1MAX,N2MAX)
          go to 41
        ENDIF
      ENDIF
 
C   ---------- Cas d'un MODELE de RAYONNEMENT
      IF (IRAYE.EQ.1) then
        if(laprop.eq.1) icavit=1
        if(laprop.eq.9) then
          call lirent( nbdang,1,iretou)
          if(ierr.ne.0) return
        endif
        if(laprop.eq.8) then
          call lirent( nbga,1,iretou)
          if(ierr.ne.0) return
        endif
        if(laprop.eq.7)then
          isyme=1
          call lirobj('POINT',ipp1,1,iretou)
          if(ierr.ne.0) return
          call lirobj('POINT',ipp2,1,iretou)
          if(idim.eq.3)call lirobj('POINT',ipp3,1,iretou)
          if(ierr.ne.0) return
        endif
        if(laprop.eq.2) then
          ifacaf=1
          call lirobj('MAILLAGE',ipfac1,1,iretou)
          if(ierr.ne.0) return
          call lirobj('MAILLAGE',ipfac2,1,iretou)
          if(ierr.ne.0) return
          call lirobj('MAILLAGE',ipfac3,1,iretou)
          if(ierr.ne.0) return
          call lirobj('MMODEL'  ,imoco ,1,iretou)
          if(ierr.ne.0) return
          call actobj('MAILLAGE',ipfac1,1)
          call actobj('MAILLAGE',ipfac2,1)
          call actobj('MAILLAGE',ipfac3,1)
          call actobj('MMODEL'  ,imoco,1)
        endif
        if (ierr.ne.0) return
      ENDIF
 
C   ---------- MODELE de SOURCE
      IF (ISRCE.EQ.1) then
C       Par DEFAUT, formulation generale (initialement "UNIFORME")
C       IF (laprop.eq.0) THEN
C         nmat=nmat+1
C         lespro(nmat)=moprop(1)
C       ENDIF
        ISRCE=ISRCE+1
      ELSEIF (ISRCE.EQ.2) THEN
        IF (lespro(nmat).EQ.'GAUSSIENNE') THEN
          IF (IDIM.EQ.1) THEN
            INTERR(1) = IDIM
            CALL ERREUR(1104)
            RETURN
          ENDIF
C         Source Gaussienne : par DEFAUT, SPHERIQUE
          IF (laprop.eq.0) THEN
            nmat=nmat+1
            lespro(nmat)=moprop(2)
          ELSE IF (laprop.eq.4) THEN
             IF (IDIM.LT.3) THEN
               INTERR(1)=IDIM
               CALL ERREUR(709)
               RETURN
             ENDIF
          ENDIF
        ENDIF
        ISRCE=ISRCE+1
      ENDIF
 
C   ---------- Cas d'un MODELE de MELANGE
      IF (LESFOR(1).EQ.'MELANGE') THEN
        CALL LIROBJ('MMODEL',ipmod,0,iOK)
        IF (IERR.NE.0) RETURN
C  ----- le melange par defaut est 'PARALLELE'
        IF (iOK.EQ.1) THEN
          CALL ACTOBJ('MMODEL',ipmod,1)
          IF (IERR.NE.0) RETURN
          IF (LAPROP.EQ.0) LAPROP=3
          LESMOD(1)=ipmod
        ENDIF
      ENDIF
C
      IF (LAPROP .EQ. 0) GOTO 42
      NMAT=NMAT+1
      LESPRO(NMAT)=MOPROP(LAPROP)
      GOTO 41
 
  42  CONTINUE
 
      IF (NMAT .NE. 0) THEN
C       on teste tout de suite l'existence de la donnee de la derivee
C       il ne faut pas de modele de materiau commencant par deri
        nmit=nmat
        do i=1,nmit
          if( lespro(i)(1:4).eq.'EPSI') then
            call erreur(19)
            return
          endif
        enddo
 
C       on cherche le mot 'EPSI'
        CALL LIRMOT(deriv,1,itrou,0)
        IF(itrou.ne.0) THEN
c           call lirmot(mderiv,5,iret,1)
c           if(ierr.ne.0) return
c           Jderiv=iret
cbp,2020-12-10 : abandon de MEPSIL (CCOPTIO) et IDERIV (MMODEL)
          MOTERR(1:40)='MODE ... EPSI ... ;'
          CALL ERREUR(1056)
          RETURN
        ENDIF
 
        IF (LESFOR(1).EQ.'THERMIQUE'.AND.ISRCE.EQ.0) THEN
C         +---------------------------------------------------------+
C         | FORMULATION THERMIQUE : 'ISOTROPE'                      |
C         +---------------------------------------------------------+
          IPROP = 3
          IF (IDIM.EQ.1) IPROP = 1
          CALL PLACE(MOPROP,IPROP,IPLAC,LESPRO(1))
          IF (IPLAC.EQ.0) THEN
            DO i=NMAT,1,-1
              LESPRO(i+1)=LESPRO(i)
            ENDDO
            LESPRO(1)='ISOTROPE'
            NMAT=NMAT+1
          ELSEif(NMAT.EQ.1)THEN
            NMAT=NMAT+1
            LESPRO(2)='CONDUCTION'
          ENDIF
 
C         Ajout du mot 'CONDUCTION' si besoin avec phase et advection
          idoico=0
          idejco=0
          DO i=1,nmat
            CALL PLACE (OPTEMP,3,iplac,LESPRO(i))
            if(iplac.eq.1.or.iplac.eq.2) idoico=1
            if(iplac.eq.3) idejco=1
          enddo
 
          if( idoico.ne.0.and.idejco.eq.0) then
            nmat=nmat+1
            lespro(nmat)='CONDUCTION'
          endif
 
        ELSEIF (LESFOR(1).EQ.'MECANIQUE'.OR. LESFOR(1).EQ.'POREUX') THEN
C         +----------------------------------------------------------+
C         | FORMULATION MECANIQUE / POREUX : 'ELASTIQUE' 'ISOTROPE'  |
C         +----------------------------------------------------------+
          IF (NMAT.GE.2)THEN
            CALL MODELA(MOPROP,NMOD)
            CALL PLACE(MOPROP,NMOD,IPLAC,LESPRO(2))
            IF (IPLAC.EQ.0) THEN
              DO i=NMAT,2,-1
                LESPRO(i+1)=LESPRO(i)
              ENDDO
              LESPRO(2)='ISOTROPE'
              NMAT=NMAT+1
            ENDIF
          ELSE IF (NMAT.EQ.1) THEN
            LESPRO(2)='ISOTROPE'
            NMAT=2
          ENDIF
C  MECANIQUE / POREUX : modele par defaut en comportement non lineaire
          CALL MODNLI(MOPROP,NMOD)
          CALL PLACE(MOPROP,NMOD,IPLAC,LESPRO(NMAT))
C  Par defaut : PLASTIQUE ISOTROPE
          IF (IPLAC.EQ.1) THEN
            NMAT=NMAT+1
            LESPRO(NMAT)='ISOTROPE'
C  Par defaut : FLUAGE NORTON
          ELSE IF (IPLAC.EQ.2) THEN
            NMAT=NMAT+1
            LESPRO(NMAT)='NORTON'
C  Par defaut : VISCOPLASTIQUE ONERA
          ELSE IF (IPLAC.EQ.3) THEN
            NMAT=NMAT+1
            LESPRO(NMAT)='ONERA'
C  Par defaut : ENDOMMAGEMENT MAZARS
          ELSE IF (IPLAC.EQ.4) THEN
            NMAT=NMAT+1
            LESPRO(NMAT)='MAZARS'
C  Par defaut : ENDOMMAGEMENT PLASTIQUE P/Y
          ELSE IF (IPLAC.EQ.5) THEN
            NMAT=NMAT+1
            LESPRO(NMAT)='PSURY'
          ELSE IF (IPLAC.EQ.6) THEN
C  Si 'MECANIQUE' OU 'POREUX' : pas de comportement par defaut
C pour 'NON_LINEAIRE'
              CALL ERREUR(945)
              RETURN
          ELSE IF (IPLAC.EQ.7) THEN
C  Si 'MECANIQUE' : pas de comportement par defaut pour 'VISCO_EXTERNE'
            IF (LESFOR(1).EQ.'MECANIQUE') THEN
              CALL ERREUR(946)
C  Si 'POREUX' : option non implementee
            ELSE IF (LESFOR(1).EQ.'POREUX') THEN
              CALL ERREUR(251)
            ENDIF
            RETURN
          ENDIF
 
        ELSEIF (LESFOR(1).EQ.'MAGNETODYNAMIQUE') THEN
C         +---------------------------------------------------------------+
C         | FORMULATION MAGNETODYNAMIQUE : 'POTENTIEL_VECTEUR' 'ISOTROPE' |
C         +---------------------------------------------------------------+
          IF (NMAT.EQ.1) THEN
            IF (LESPRO(1).NE.'POTENTIEL_VECTEU') THEN
              LESPRO(2)=LESPRO(1)
              LESPRO(1)='POTENTIEL_VECTEU'
            ELSE
              LESPRO(2)='ISOTROPE'
            ENDIF
            NMAT=2
          ENDIF
 
        ELSEIF (LESFOR(1).EQ.'MELANGE'         ) THEN
C         +-------------------------------+
C         | FORMULATION MELANGE : 'CEREM' |
C         +-------------------------------+
          NMAT1=NMAT
          IF (NMAT.EQ.0) THEN
            LESPRO(1)='CEREM'
            NMAT=1
          ENDIF
 
        ELSEIF (LESFOR(1).EQ.'LIAISON'         ) THEN
C         +-------------------------------------------------+
C         | FORMULATION  LIAISON : pas d''option par defaut |
C         +-------------------------------------------------+
 
        ELSEIF (LESFOR(1).EQ.'ELECTROSTATIQUE' ) THEN
C         +-------------------------------------------+
C         | FORMULATION  ELECTROSTATIQUE : 'ISOTROPE' |
C         +-------------------------------------------+
          IPROP = 3
          IF (IDIM.EQ.1) IPROP = 1
          CALL PLACE(MOPROP(1),IPROP,IPLAC,LESPRO(1))
          IF (IPLAC.EQ.0) THEN
            DO i=NMAT,1,-1
              LESPRO(i+1)=LESPRO(i)
            ENDDO
            LESPRO(1)='ISOTROPE'
            NMAT=NMAT+1
          ENDIF
 
        ELSEIF (LESFOR(1).EQ.'DIFFUSION'       ) THEN
C         +-------------------------------------------+
C         | FORMULATION DIFFUSION : 'ISOTROPE' 'FICK' |
C         +-------------------------------------------+
          IPROP = 3
          IF (IDIM.EQ.1) IPROP = 1
          CALL PLACE(MOPROP(1),IPROP,IPLAC,LESPRO(1))
          IF (IPLAC.EQ.0) THEN
            DO i=NMAT,1,-1
              LESPRO(i+1)=LESPRO(i)
            ENDDO
            LESPRO(1)='ISOTROPE'
            NMAT=NMAT+1
          ENDIF
          CALL MODDIF(MOPROP,NMOD)
          CALL PLACE(MOPROP,NMOD,IPLAC,LESPRO(NMAT))
          IF (IPLAC.EQ.0) THEN
            NMAT=NMAT+1
            LESPRO(NMAT)='FICK'
          ENDIF
 
C         Ajout du mot 'FICK' si besoin avec 'ADVECTION'
          CALL PLACE(LESPRO,nmat,iplac,'ADVECTION')
          if(iplac .gt. 0) then
            NMAT=NMAT+1
            LESPRO(NMAT)='FICK'
          endif
 
        ELSEIF (LESFOR(1).EQ.'CONTACT'         ) THEN
C         +----------------------------------+
C         | FORMULATION CONTACT : UNILATERAL |
C         +----------------------------------+
          call place ( moprop,2,iplac,lespro(1))
          if( iplac.eq.0) then
           do iur=1,nmat
             lespro(nmat+2-iur)=lespro (nmat +1-iur)
           enddo
           lespro(1)='UNILATERAL'
           nmat=nmat+1
          endif
 
        ELSEIF (LESFOR(1).EQ.'CONTRAINTE'      ) THEN
C         +----------------------------------+
C         | FORMULATION CONTRAINTE           |
C         +----------------------------------+
          call place ( moprop,3,iplac,lespro(1))
          if (iplac.eq.0) then
           do iur=1,nmat
             lespro(nmat+2-iur)=lespro (nmat +1-iur)
           enddo
           lespro(1)='CINEMATIQUE'
           nmat=nmat+1
          endif
        ENDIF
 
      ELSE
C       si NMAT=0 on met le premier mot autorise
        NMAT     = 1
        LESPRO(1)= MOPROP(1)
 
        IF    (LESFOR(1).EQ.'CHARGEMENT') THEN
C         +------------------------------------------------------------------+
C         | Defaut pour une FORMULATION CHARGEMENT : PAS DE CHOIX PAR DEFAUT |
C         +------------------------------------------------------------------+
C         L'UTILISATEUR DOIT SPECIFIER D'AUTRES MOT CLES APRES 'CHARGEMENT'
          CALL ERREUR(251)
          RETURN
 
        ELSEIF(LESFOR(1).EQ.'THERMIQUE') THEN
C         +----------------------------------------------------+
C         | Defaut pour une FORMULATION THERMIQUE : CONDUCTION |
C         +----------------------------------------------------+
          NMAT  = NMAT+1
          LESPRO(NMAT)='CONDUCTION'
 
        ELSEIF (LESFOR(1).EQ.'MECANIQUE'.OR.
     &          LESFOR(1).EQ.'POREUX'   .OR.
     &          LESFOR(1).EQ.'MAGNETODYNAMIQUE') THEN
C         +------------------------------------------------------------------------------+
C         | Defaut pour une FORMULATION MECANIQUE, POREUX ou MAGNETODYNAMIQUE : ISOTROPE |
C         +------------------------------------------------------------------------------+
          NMAT=NMAT+1
          LESPRO(NMAT)='ISOTROPE'
 
        ELSEIF (LESFOR(1).EQ.'DIFFUSION') THEN
C         +----------------------------------------------+
C         | Defaut pour une FORMULATION DIFFUSION : FICK |
C         +----------------------------------------------+
          NMAT  = NMAT+1
          LESPRO(NMAT)='FICK'
 
        ELSEIF (LESFOR(1).EQ.'NAVIER_STOKES'.OR.
     &           LESFOR(1).EQ.'EULER') THEN
C         +----------------------------------------------------------------+
C         | Defaut pour une FORMULATION NAVIER_STOKES OU EULER : NEWTONIEN |
C         +----------------------------------------------------------------+
          NMAT = 1
          LESPRO(NMAT)='NEWTONIEN'
 
        ELSEIF (LESFOR(1).EQ.'FISSURE') THEN
C         +-------------------------------------+
C         | Defaut pour une FORMULATION FISSURE |
C         +-------------------------------------+
          NMAT = 3
          LESPRO(1)='MASS'
          LESPRO(2)='PARF'
          LESPRO(3)='POISEU_BLASIUS'
 
        ELSEIF(LESFOR(1).EQ.'CONTACT') THEN
C         +---------------------------------------------------+
C         | Defaut pour une FORMULATION CONTACT : UNILATERAL  |
C         +---------------------------------------------------+
          NMAT=1
          LESPRO(1)='UNILATERAL'
 
        ELSEIF(LESFOR(1).EQ.'CONTRAINTE') THEN
C         +---------------------------------------------------+
C         | Defaut pour une FORMULATION CONTRAINTE            |
C         +---------------------------------------------------+
          NMAT=1
          LESPRO(1)='CINEMATIQUE'
        ENDIF
      ENDIF
 
      IF(LESFOR(1).EQ.'CHANGEMENT_PHASE' ) THEN
C       +------------------------------------------------------------------------+
C       | FORMULATION CHANGEMENT_PHASE : LECTURE DES INCONNUES PRIMALES & DUALES |
C       +------------------------------------------------------------------------+
        CALL LIRMOT(MOINCO,NBDIF,IPLAC,0)
        IF (IPLAC.EQ.0) THEN
          CALL ERREUR(1093)
          RETURN
 
        ELSE
          IF     (LESPRO(1)(1:10).EQ.'PARFAIT   ') THEN
            JGM=2
          ELSEIF (LESPRO(1)(1:10).EQ.'SOLUBILITE') THEN
            JGM=4
          ELSE
            CALL ERREUR(5)
          ENDIF
          JGN   =LOCOMP
          SEGINI,MLMOT1
          IPRIDU=MLMOT1
          DO IMOT=1,JGM
            CALL LIRCHA(MOPRID,1,ILONG)
            IF (IERR.NE.0) RETURN
            MLMOT1.MOTS(IMOT) = MOPRID
          ENDDO
        ENDIF
      ENDIF
 
 43   CONTINUE
C  Lecture eventuelle des types d'ELEMENTS FINIS a utiliser
      ITEF=0
C=DEB==== FORMULATION HHO ==== Cas particulier =========================
      CALL LIRMOT(mcHHO,NMHHO,iHHO,0)
      IF (IERR.NE.0) RETURN
      IF (iHHO.NE.0) THEN
        CALL REFUS
        CALL LIRCHA(chaHHO,1,IRETI)
        IF (IERR.NE.0) RETURN
        loHHO = .TRUE.
      END IF
C=FIN==== FORMULATION HHO ==============================================
      IF (NBTEF.EQ.0) GOTO 2
C      WRITE(*,*) 'MODELI:',(MOTEF(i),':',i=1,NBTEF)
      CALL MESLIR(-178)
 1    CONTINUE
      CALL LIRMOT(MOTEF,NBTEF,LETEF,0)
      IF (IERR.NE.0) RETURN
      IF (LETEF.EQ.0) GOTO 2
      ITEF=ITEF+1
      LESTEF(ITEF)=MOTEF(LETEF)
      if (ITEF.eq.1.and.lesfor(1).eq.'NAVIER_STOKES') goto 2
      CALL MESLIR(-177)
      GOTO 1
 
 2    CONTINUE
c  Lecture eventuelle de listmots
      jlmot1 = 0
      jlmot2 = 0
      CALL LIROBJ('LISTMOTS',jlmot1,0,iret)
      if (ierr.ne.0) return
      if (jlmot1.gt.0) then
        call lirobj('LISTMOTS',jlmot2,1,iret)
        if (ierr.ne.0) return
        mlmot5 = jlmot1
        mlmot6 = jlmot2
        segact,mlmot5,mlmot6
        if (mlmot5.mots(/2).ne.mlmot6.mots(/2)) then
          call erreur(26)
          return
        endif
      endif
 
C   lecture pour mecanique aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa
C  En formulations 'MECANIQUE' et 'POREUX' : detection d'une loi non
C  lineaire externe, le cas echeant saisie de donnees complementaires.
C  Caracterisation : loi non lineaire externe si
C  - famille 'VISCO_EXTERNE' ou
C  - famille 'NON_LINEAIRE', materiau 'UTILISATEUR'.
C  si pas loi externe lecture eventuelle des parametres externes
      LMEEXT=.FALSE.
      LMEVIX=.FALSE.
      LMENLX=.FALSE.
      LMENUM = 0
      LMENOM = ' '
      LMELIB = ' '
      LMEFCT = ' '
      LMELGB = 0
      LMELGT = 0
      LMELOI = 0
      LMEPTR = 0
 
      IF ( (NFOR.EQ.1).AND.
     &  (LESFOR(1).EQ.'MECANIQUE'.OR.LESFOR(1).EQ.'POREUX') ) THEN
        DO i=1,NMAT
          IF (LESPRO(i).EQ.'VISCO_EXTERNE') THEN
            LMEVIX=.TRUE.
            GOTO 203
          ENDIF
        ENDDO
 
        IF (.NOT.LMEVIX) THEN
          DO i=1,NMAT
            IF (LESPRO(i).EQ.'UTILISATEUR') THEN
              LMENLX = .TRUE.
              GOTO 203
            ENDIF
          ENDDO
        ENDIF
C........N.B. LMEEXT exprime la condition (NFOR.EQ.1) ET
C             (LESFOR(1).EQ.'MECANIQUE') ET (loi non lineaire externe)
 203     LMEEXT = LMEVIX.OR.LMENLX
         IF ( LMEEXT ) THEN
C  lecture et verif des noms des materiaux, des
C  noms des variables internes,  des noms des parametre externe pour
C loi externes
 210        CALL LIRMOT(MOEXT,NBEXT,LEXT,0)
C           Si on ne trouve plus l'un des mots cles attendus, on sort
            IF (LEXT.EQ.0) GOTO 211
C           Lecture d'un entier sous 'NUME_LOI'
            IF (LEXT.EQ.1) THEN
               CALL LIRENT(LMENUM,1,IRET)
               IF (IERR.NE.0) RETURN
C           Valeur illicite du numero de la loi (superieur ou egal a 1)
               IF (LMENUM.LT.1 .OR. LMENUM.GE.1000000) THEN
                 INTERR(1) = LMENUM
                 CALL ERREUR(36)
                 CALL ERREUR(947)
                 RETURN
               ENDIF
C           Lecture du nom de la loi sous 'NOM_LOI'
            ELSE IF (LEXT.EQ.2) THEN
               MOTEMP = ' '
               CALL LIRCHA(MOTEMP,1,IRET)
               IF (IERR.NE.0) RETURN
               IRET = LONG(MOTEMP(1:IRET))
               IF (IRET.GT.16) THEN
                 INTERR(1) = IRET
                 MOTERR = MOTEMP(1:IRET)
                 CALL ERREUR(-2)
                 CALL ERREUR(21)
                 RETURN
               ELSE IF (IRET.LE.0) THEN
                 INTERR(1) = 0
                 MOTERR = 'NOM_LOI'
                 CALL ERREUR(-2)
                 CALL ERREUR(6)
                 RETURN
               ENDIF
               LMENOM = '        '
               LMENOM(1:IRET) = MOTEMP(1:IRET)
C           Lecture d'un objet LISTMOTS sous 'PARA_LOI'
            ELSE IF (LEXT.EQ.3) THEN
               CALL LIROBJ('LISTMOTS',LUPARX,1,IRET)
               IF (IERR.NE.0) RETURN
C           Lecture d'un objet LISTMOTS sous 'C_MATERIAU'
            ELSE IF (LEXT.EQ.4) THEN
               CALL LIROBJ('LISTMOTS',LUCMAT,1,IRET)
               IF (IERR.NE.0) RETURN
C           Lecture d'un objet LISTMOTS sous 'C_VARINTER'
            ELSE IF (LEXT.EQ.5) THEN
               CALL LIROBJ('LISTMOTS',LUCVAR,1,IRET)
               IF (IERR.NE.0) RETURN
C           Lecture du nom (du fichier) de la bibliotheque de la loi
            ELSE IF (LEXT.EQ.6) THEN
               MOTEMP = ' '
               CALL LIRCHA(MOTEMP,1,IRET)
               IF (IERR.NE.0) RETURN
               IF (IRET.GT.LOCHAI) THEN
                 CALL ERREUR(1110)
                 RETURN
               ENDIF
               IRET = LONG(MOTEMP(1:IRET))
               IF (IRET.LE.0) THEN
                 INTERR(1) = 0
                 MOTERR = 'LIB_LOI'
                 CALL ERREUR(-2)
                 CALL ERREUR(6)
                 RETURN
               END IF
               LMELIB = '        '
               LMELIB(1:IRET) = MOTEMP(1:IRET)
               LMELGB = IRET
               LMEPTR = IRET
C           Lecture du nom de la fonction de la loi
            ELSE IF (LEXT.EQ.7) THEN
               MOTEMP = ' '
               CALL LIRCHA(MOTEMP,1,IRET)
               IF (IERR.NE.0) RETURN
               IF (IRET.GT.LOCHAI) THEN
                 CALL ERREUR(1110)
                 RETURN
               ENDIF
               IRET = LONG(MOTEMP(1:IRET))
               IF (IRET.LE.0) THEN
                 INTERR(1) = 0
                 MOTERR = 'FCT_LOI'
                 CALL ERREUR(-2)
                 CALL ERREUR(6)
                 RETURN
               ENDIF
               LMEFCT = ' '
               LMEFCT(1:IRET) = MOTEMP(1:IRET)
               LMELGT = IRET
            ENDIF
C           On repete jusqu'a ce qu'on ne trouve plus aucun des
C           mots cles attendus, regle de surcharge le cas echeant
            GOTO 210
 211        CONTINUE
C...........Verifications sur les donnees
C           Il manque 'NUME_LOI' ou 'NOM_LOI' (toujours obligatoire)
            IF (LMENUM.EQ.0 .AND. LMENOM.EQ.'                ') THEN
              IF (LMELGT.EQ.0) THEN
                CALL ERREUR(641)
                RETURN
              ENDIF
            ENDIF
            IF (LMENUM.NE.0 .AND. LMENOM.NE.'                ') THEN
               MOTERR(1:16) = 'NUME_LOINOM_LOI '
               CALL ERREUR(135)
               RETURN
            ENDIF
C Les liste des composantes ne doivent pas etre vides.
            DO i = 1, 3
              IF (i.EQ.1) mlmots = LUPARX
              IF (i.EQ.2) mlmots = LUCMAT
              IF (i.EQ.3) mlmots = LUCVAR
              IF (mlmots.NE.0) THEN
                SEGACT,mlmots
                NBCOMP = mlmots.mots(/2)
                IF (NBCOMP.EQ.0) THEN
                  CALL ERREUR(964)
                  RETURN
                ENDIF
              ENDIF
            ENDDO
C Dans le cas d'un modele NON_LINEAIRE UTILISATEUR, on rajoute en fin de
C liste des proprietes du modele, le numero ou le nom de la loi attribue
C par l'utilisateur.
            NMAT = NMAT + 1
            LESPRO(NMAT) = '                '
            IF (LMENUM.EQ.0) THEN
              LESPRO(NMAT) = LMENOM
              IF (LMELGT.GT.0 .AND. LMENOM.EQ.'                ') THEN
c* On espere mettre un numero "unique" dans le nom !
                SEGINI,limode
                WRITE(LESPRO(NMAT)(1:16),'(I16)') limode
                SEGSUP,limode
              ENDIF
            ELSE
              WRITE(LESPRO(NMAT)(1:16),'(I16)') LMENUM
            ENDIF
C           Verifications pour une loi 'NON_LINEAIRE' 'UTILISATEUR'
            IF ( LMENLX ) THEN
C              Il manque les composantes materielles sous 'C_MATERIAU'
               IF (LUCMAT.EQ.0) THEN
                  CALL ERREUR(641)
                  RETURN
               ENDIF
C              La liste des composantes materielles saisie sous
C              'C_MATERIAU' ne doit pas etre vide
               MLMOTS=LUCMAT
               SEGACT,MLMOTS
               NBCOMP = MOTS(/2)
               IF (NBCOMP.EQ.0) THEN
                  CALL ERREUR(964)
                  RETURN
               ENDIF
            ENDIF
C Dans le cas d'une libraire externe, quelques verifications puis
C recherche du pointeur de la fonction externe
            IF (LMEPTR.GT.0) THEN
C Si le nom de la fonction n'a pas ete fourni avec le mot-cle 'FCT_LOI',
C on le construit a partir de 'NOM_LOI' ou 'NUME_LOI'.
              IF (LMELGT.EQ.0) THEN
                LMEFCT = ' '
                IF (LMENUM.EQ.0) THEN
                  IRET = LONG(LMENOM)
                  LMEFCT(1:IRET) = LMENOM(1:IRET)
                  LMELGT = IRET
                ELSE
                  IRET = 0
                  DO i = 1, 16
                    IRET = IRET + 1
                    IF (LESPRO(NMAT)(i:i).NE.' ') GOTO 220
                  ENDDO
 220              continue
                  LMEFCT = 'umat_'//LESPRO(NMAT)(IRET:16)
                  LMELGT = 22-IRET
                ENDIF
              ENDIF
              ip = -1
              CALL LEXTOP(LMELIB,LMEFCT,ip,LMELOI,LMEPTR)
***           IF (IERR.NE.0) RETURN
*si pas d'erreur LMELOI > 0 et LMEPTR >0 pointe sur une fonction
*dbg              IF (LMELOI.LE.0) CALL ERREUR(5)
*dbg              IF (LMEPTR.LE.0) CALL ERREUR(5)
              LMELGB = LONG(LMELIB)
              LMELGT = LONG(LMEFCT)
*dbg            write(ioimp,*) 'LMELOI =',LMELOI,LMEPTR,LMELGB,LMELGT,
*dbg     &                     LMELIB(1:LMELGB),'=',LMEFCT(1:LMELGT)
            ENDIF
         ELSE
C  si pas lois externes lecture des noms des parametres externes
           CALL LIRMOT(MOEXT(2),1,LEXT,0)
           IF (LEXT.NE.0) THEN
             CALL LIROBJ('LISTMOTS',luparx,1,iret)
             IF (IERR.NE.0) RETURN
           ENDIF
         ENDIF
C           Verifications sur les parametres, si declares
         IF (luparx.GT.0) THEN
C     Si la temperature 'T   ' fait partie des parametres de
C     la loi, elle doit etre declaree en tete
            mlmots=luparx
            SEGACT,MLMOTS
            NBPARA=MOTS(/2)
            IF (NBPARA.GT.0) THEN
              DO IP = 1, NBPARA
                IF (MOTS(IP).EQ.'T   ') THEN
                  IF (IP.GT.1) THEN
                    CALL ERREUR(948)
                    RETURN
                  ENDIF
                  GOTO 221
                ENDIF
              ENDDO
 221          CONTINUE
            ENDIF
C     Pas de parametres redondants
            IF (NBPARA.GT.1) THEN
               DO 230 IP1 = 1, NBPARA-1
                  PAR1 = MOTS(IP1)
                  DO 231 IP2 = IP1+1, NBPARA
                     IF (MOTS(IP2).EQ.PAR1) THEN
                        CALL ERREUR(949)
                        RETURN
                     ENDIF
 231              CONTINUE
 230           CONTINUE
            ENDIF
         ENDIF
      ENDIF
C  fin lecture mecanique aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa
 
C Formulation 'DIFFUSION' : oooooooooooooooooooooooooooooooooooooooooooo
C  - Detection d'une loi non lineaire externe (mot-cle 'UTILISATEUR')
C  - Detection de la quantite de l'effet Soret (mot-cle 'SORET')
C  - Lecture de l'inconnue "diffusant" (mot-cle 'INCO')
C 3    CONTINUE
      LDIEXT = .FALSE.
      LDISOR = .FALSE.
      LDINUM = 0
      LDINOM = ' '
      LDILIB = ' '
      LDIFCT = ' '
      LDILGB = 0
      LDILGT = 0
      LDILOI = 0
      LDIPTR = 0
      IF (NFOR.EQ.1 .AND. LESFOR(1).EQ.'DIFFUSION') THEN
C  -- Recherche des informations sur la presence d'une loi externe --
        DO i=1,NMAT
          LDIEXT = LESPRO(i).EQ.'UTILISATEUR'
          LDISOR = LESPRO(i).EQ.'SORET'
        ENDDO
C - Lecture des informations pour la loi externe
        IF (LDIEXT) THEN
 310      CONTINUE
          CALL LIRMOT(MOEXT,NBEXT,LEXT,0)
          IF (LEXT.EQ.0) GOTO 311
C       Lecture d'un entier sous 'NUME_LOI'
          IF (LEXT.EQ.1) THEN
            CALL LIRENT(LDINUM,1,IRET)
            IF (IERR.NE.0) RETURN
            IF (LDINUM.LT.1 .OR. LDINUM.GE.1000000) THEN
              INTERR(1) = LDINUM
              CALL ERREUR(36)
              CALL ERREUR(947)
              RETURN
            ENDIF
C           Lecture du nom de la loi sous 'NOM_LOI'
          ELSE IF (LEXT.EQ.2) THEN
            MOTEMP = ' '
            CALL LIRCHA(MOTEMP,1,IRET)
            IF (IERR.NE.0) RETURN
            IRET = LONG(MOTEMP(1:IRET))
            IF (IRET.GT.16) THEN
              INTERR(1) = IRET
              MOTERR = MOTEMP(1:IRET)
              CALL ERREUR(-2)
              CALL ERREUR(21)
              RETURN
            ELSE IF (IRET.LE.0) THEN
              INTERR(1) = IRET
              MOTERR = 'NOM_LOI'
              CALL ERREUR(-2)
              CALL ERREUR(6)
              RETURN
            ENDIF
            LDINOM = ' '
            LDINOM(1:IRET) = MOTEMP(1:IRET)
C        Lecture d'un objet LISTMOTS sous 'PARA_LOI'
          ELSE IF (LEXT.EQ.3) THEN
            CALL LIROBJ('LISTMOTS',luparx,1,IRET)
            IF (IERR.NE.0) RETURN
C       Lecture d'un objet LISTMOTS sous 'C_MATERIAU'
          ELSE IF (LEXT.EQ.4) THEN
            CALL LIROBJ('LISTMOTS',lucmat,1,IRET)
            IF (IERR.NE.0) RETURN
C       Lecture d'un objet LISTMOTS sous 'C_VARINTER'
          ELSE IF (LEXT.EQ.5) THEN
            CALL LIROBJ('LISTMOTS',lucvar,1,IRET)
            IF (IERR.NE.0) RETURN
C       Lecture du nom (du fichier) de la bibliotheque de la loi
          ELSE IF (LEXT.EQ.6) THEN
             MOTEMP = ' '
             CALL LIRCHA(MOTEMP,1,IRET)
             IF (IERR.NE.0) RETURN
             IF (IRET.GT.LOCHAI) THEN
               CALL ERREUR(1110)
               RETURN
             ENDIF
             IRET = LONG(MOTEMP(1:IRET))
             IF (IRET.LE.0) THEN
               INTERR(1) = IRET
               MOTERR = MOTEMP
               CALL ERREUR(-2)
               CALL ERREUR(6)
               RETURN
             ENDIF
             LDILIB = ' '
             LDILIB(1:IRET) = MOTEMP(1:IRET)
             LDILGB = IRET
             LDIPTR = IRET
C           Lecture du nom de la fonction de la loi
          ELSE IF (LEXT.EQ.7) THEN
             MOTEMP = ' '
             CALL LIRCHA(MOTEMP,1,IRET)
             IF (IERR.NE.0) RETURN
             IF (IRET.GT.LOCHAI) THEN
               CALL ERREUR(1110)
               RETURN
             ENDIF
             IRET = LONG(MOTEMP(1:IRET))
             IF (IRET.LE.0) THEN
               INTERR(1) = IRET
               MOTERR = MOTEMP
               CALL ERREUR(-2)
               CALL ERREUR(6)
               RETURN
             ENDIF
             LDIFCT(1:IRET) = MOTEMP(1:IRET)
             LDILGT = IRET
          ENDIF
          GOTO 310
 311      CONTINUE
C Verifications des informations obligatoires de la loi externe
C       Il manque 'NUME_LOI' ou 'NOM_LOI' (toujours obligatoire)
          IF (LDINUM.EQ.0 .AND. LDINOM.EQ.'                ') THEN
            if (LDILGT.eq.0) then
              CALL ERREUR(641)
              RETURN
            endif
          ENDIF
          IF (LDINUM.NE.0 .AND. LDINOM.NE.'                ') THEN
            MOTERR(1:16) = 'NUME_LOINOM_LOI '
            CALL ERREUR(135)
            RETURN
          ENDIF
C       Il manque la liste 'C_MATERIAU'
          IF (lucmat.EQ.0) THEN
            CALL ERREUR(641)
            RETURN
          ENDIF
C Les liste des composantes ne doivent pas etre vides.
          DO i = 1, 3
            IF (i.EQ.1) MLMOTS = luparx
            IF (i.EQ.2) MLMOTS = lucmat
            IF (i.EQ.3) MLMOTS = lucvar
            IF (MLMOTS.NE.0) THEN
              SEGACT,MLMOTS
              NBCOMP = MOTS(/2)
              IF (NBCOMP.EQ.0) THEN
                CALL ERREUR(964)
                RETURN
              ENDIF
            ENDIF
          ENDDO
C Dans le cas d'un modele UTILISATEUR, on rajoute en fin de
C liste des proprietes du modele, le numero attribue par l'utilisateur.
          NMAT = NMAT + 1
          LESPRO(NMAT) = '                '
          IF (LDINUM.EQ.0) THEN
            LESPRO(NMAT) = LDINOM
            if (LDILGT.gt.0.and.LDINOM.eq.'                ') then
c* On espere mettre un numero "unique" dans le nom !
              segini,limode
              write(LESPRO(NMAT)(1:16),'(I16)') limode
              segsup,limode
            endif
          ELSE
            WRITE(LESPRO(NMAT)(1:16),'(I16)') LDINUM
          ENDIF
C Dans le cas d'une libraire externe, quelques verifications puis
C recherche du pointeur de la fonction externe
          IF (LDIPTR.GT.0) THEN
C Si le nom de la fonction n'a pas ete fourni avec le mot-cle 'FCT_LOI',
C on le construit a partir de 'NOM_LOI' ou 'NUME_LOI'.
           IF (LDILGT.EQ.0) THEN
            LDIFCT = ' '
            IF (LDINUM.EQ.0) THEN
              IRET = LONG(LDINOM)
              LDIFCT(1:IRET) = LDINOM(1:IRET)
              LDILGT = IRET
            ELSE
              IRET = 0
              DO i = 1, 16
                IRET = IRET + 1
                IF (LESPRO(NMAT)(i:i).NE.' ') GOTO 320
              ENDDO
 320          CONTINUE
              LDIFCT = 'umat_'//LESPRO(NMAT)(IRET:16)
              LDILGT = 22-IRET
            ENDIF
           ENDIF
            ip = -1
            CALL LEXTOP(LDILIB,LDIFCT,ip,LDILOI,LDIPTR)
            IF (IERR.NE.0) RETURN
*si pas d'erreur LDILOI > 0 et LDIPTR >0 pointe sur une fonction
*dbg            IF (LDILOI.LE.0) CALL ERREUR(5)
*dbg            IF (LDIPTR.LE.0) CALL ERREUR(5)
            LDILGB = LONG(LDILIB)
            LDILGT = LONG(LDIFCT)
          ENDIF
        ENDIF
C - Lecture des informations pour la loi Soret :
C - quantite dont le gradient est l'origine de l'effet ('T' par defaut)
        IF (LDISOR) THEN
          mlmots = 0
          CHARIN = 'T       '
C   Lecture du mot-cle 'PARA_LOI' et donnees associees
          CALL LIRMOT(MOEXT(2),1,LEXT,0)
          IF (IERR.NE.0) RETURN
          IF (LEXT.EQ.1) THEN
            CALL LIROBJ('LISTMOTS',mlmots,0,IRET)
            IF (IERR.NE.0) RETURN
            IF (IRET.EQ.0) THEN
              CALL LIRCHA(CHARIN,1,IRETI)
              IF (IERR.NE.0) RETURN
              IRETI=LONG(CHARIN)
              IF (IRETI.EQ.0) CALL ERREUR(643)
            ELSE
              SEGACT,mlmots
              NBCOMP = mots(/2)
              IF (NBCOMP.EQ.0) THEN
                CALL ERREUR(964)
              ELSE
                CHARIN = MOTS(1)
                IRETI = LONG(CHARIN)
                IF (IRETI.EQ.0) CALL ERREUR(643)
              ENDIF
            ENDIF
            IF (IERR.NE.0) RETURN
            IRETMA = 2
C*8         IRETMA = 6
            IF (IRETI.GT.IRETMA) THEN
              INTERR(1)   = IRETMA
              MOTERR(1:8) = CHARIN(1:IRETI)
              CALL ERREUR(-353)
            ENDIF
            IRETI = MIN(IRETI,IRETMA)
            CHARIN(IRETI+1:8) = ' '
          ENDIF
          JGM = 1
          JGN = LOCOMP
          SEGINI,mlmots
          mots(1) = CHARIN
          luparx = mlmots
        ENDIF
 
C  -- Pour la formulation DIFFUSION : lecture quantite (ddl) diffusant --
C  -- On cherche a lire le mot 'INCO' suivi du nom de l'INCOnnue donne --
C  -- soit par un LISTMOTS, soit par un MOT puis eventuellement du nom --
C  -- de la grandeur DUALe donne par un objet de meme type que pour le --
C  -- nom de l'inconnue.                                               --
        CALL LIRMOT(MOINCO,NBDIF,LEXT,0)
        IF (LEXT.EQ.0) THEN
C*8       MDIINC='CONC    '
C*8       MDIDUA='QCONC   '
          MDIINC='CO      '
          MDIDUA='QCO     '
        ELSE
          MDIINC='        '
          MDIDUA='Q       '
          CHARIN='        '
          CHARRE='        '
          CALL LIROBJ('LISTMOTS',mlmots,0,IRET)
          IF (IERR.NE.0) RETURN
          IF (mlmots.NE.0) THEN
            SEGACT,mlmots
            NBCOMP = MOTS(/2)
            IF (NBCOMP.EQ.0) THEN
              CALL ERREUR(643)
            ELSE
              CHARIN=MOTS(1)
              IRETI=LONG(CHARIN)
              IF (IRETI.EQ.0) CALL ERREUR(643)
            ENDIF
            IF (IERR.NE.0) RETURN
            CALL LIROBJ('LISTMOTS',mlmots,0,IRETE)
            IF (IERR.NE.0) RETURN
            IF (mlmots.NE.0) THEN
              SEGACT,mlmots
              NBCOMP = MOTS(/2)
              IF (NBCOMP.EQ.0) THEN
                CALL ERREUR(643)
              ELSE
                CHARRE=MOTS(1)
                IRETE=LONG(CHARRE)
                IF (IRETE.EQ.0) CALL ERREUR(643)
              ENDIF
              IF (IERR.NE.0) RETURN
            ENDIF
 
          ELSE
            CALL LIRCHA(CHARIN,1,IRETI)
            IF (IERR.NE.0) RETURN
            IRETI = LONG(CHARIN(1:IRETI))
            IF (IRETI.EQ.0) THEN
              CALL ERREUR(643)
              RETURN
            ENDIF
            CALL LIRCHA(CHARRE,0,IRETE)
            IF (IERR.NE.0) RETURN
            IF (IRETE.GT.0) THEN
              IRETE = LONG(CHARRE(1:IRETE))
              IF (IRETE.EQ.0) THEN
                CALL ERREUR(643)
                RETURN
              ENDIF
            ENDIF
          ENDIF
 
          IRETMA = 2
C*8       IRETMA = 6
          IF (IRETI.GT.IRETMA) THEN
            INTERR(1)   = IRETMA
            MOTERR(1:8) = CHARIN(1:IRETI)
            CALL ERREUR(-353)
          ENDIF
          IRETI = MIN(IRETI,IRETMA)
          MDIINC(1:IRETI)=CHARIN(1:IRETI)
          IF (IRETE.EQ.0) THEN
            MDIDUA(2:1+IRETI)=MDIINC(1:IRETI)
 
          ELSE
            IRETMA = IRETMA + 2
            IF (IRETE.GT.IRETMA) THEN
              INTERR(1)   = IRETMA
              MOTERR(1:8) = CHARRE(1:IRETE)
              CALL ERREUR(-353)
            ENDIF
            IRETE=MIN(IRETE,IRETMA)
            MDIDUA(1:IRETE)=CHARRE(1:IRETE)
          ENDIF
        ENDIF
 
c*      Verification des noms de primale et duale lues
        CALL VERMDI(MDIINC,MDIDUA)
        IF (IERR.NE.0) RETURN
 
C       On les place dans un LISTMOTS pour TYMODE et IVAMODE
        JGN = LOCOMP
        JGM = 2
        SEGINI,MLMOT1
        iplrdi=MLMOT1
        MLMOT1.MOTS(1) = MDIINC
        MLMOT1.MOTS(2) = MDIDUA
      ENDIF
C Fin Formulation 'DIFFUSION' oooooooooooooooooooooooooooooooooooooooooo
 
C  Lecture eventuelle du NOM de CONSTITUANT, du nombre de POINTs
C  d'INTEGRATION, du point support pour les modes en DEFOrmations
C  PLANEs GENEralisees, du nom de la phase, de la formulation non_locale
C  fin des lecture en 22
 674  CONTINUE
      kcons  = 0
      NGINT  = 0
      NGRIG  = 0
      NGMAS  = 0
      NGCON  = 0
      IPTGEN = 0
      IPMOD1 = 0
      klcon  = 0
      plicon = 0
      ILIE   = 0
      INLOC  = 0
      INLVIA = 0
      LULVIA = 0
      kbnlin = 0
 
 675  CONTINUE
      CALL LIRMOT(MOCON,NBCON,LECON,0)
      IF (LECON.EQ.0) GOTO 22
      IF (LECON.EQ.1) THEN
        CALL LIRCHA(CONM,1,kcons)
        IF (IERR.NE.0) RETURN
      ELSE IF (LECON.EQ.2) THEN
        i1foi = 1
 677    continue
        legaus=0
        CALL LIRMOT(MGAUSS,4,legaus,0)
        if (ierr.ne.0) return
        if (i1foi.ne.1.and.legaus.eq.0) goto 675
        CALL LIRENT(itt,1,iret)
        if (ierr.ne.0) return
        if (itt.lt.1) then
          interr(1) = itt
          call erreur(36)
          return
        endif
        if (legaus.eq.0 .or. legaus.eq.1) then
c itt doit etre impair (> 0)
          IF (MOD(itt,2).EQ.0) THEN
            call erreur(607)
            return
          ENDIF
          NGINT = itt
        endif
        if (legaus.eq.2) NGRIG = itt
        if (legaus.eq.3) NGMAS = itt
        if (legaus.eq.4) NGCON = itt
        if (i1foi.eq.1.and.legaus.eq.0) goto 675
        i1foi = 0
c INTE itt <=> INTE EPAI itt ; autres mots a ecrire
c Syntaxe de modeli non decrite :
c Si plusieurs mots de MGAUSS 
c INTE MOT1 itt1 MOT2 itt2 ... ; (couples MOTi iiti obligatoires)
        goto 677
      ELSE IF (LECON.EQ.3) THEN
        CALL LIROBJ('POINT',IPTGEN,1,IRET)
        IF (IERR.NE.0) RETURN
C On transforme le point en maillage de POI1 (avec un seul element)
        CALL CRELEM(IPTGEN)
C On verifie s'il n'a pas deja ete preconditionne.
        CALL CRECH1(IPTGEN,1)
      ELSE IF (LECON.EQ.4) THEN
        CALL LIRCHA(PHAM,1,IRET)
        IF (IERR.NE.0) RETURN
      ELSE IF (LECON.EQ.5) THEN
        NMAT = NMAT + 1
        LESPRO(NMAT) = 'STATIONNAIRE'
        if (iptabs.gt.0) then
        else
         CALL LIROBJ('MMODEL',IPMOD1,1,IRET)
         IF (IERR.NE.0) RETURN
        endif
C LCOI et LCOS : options non documentees pour le modele LIAISON !
C Lecture obligatoire du modele associe (sinon options sans interet)
      ELSE IF (LECON.EQ.6.OR.LECON.EQ.7) THEN
        IF (LESFOR(1).NE.'LIAISON') THEN
          CALL ERREUR(251)
          RETURN
        ENDIF
        CALL LIROBJ('MMODEL  ',ipmod2,1,iret)
        IF (IERR.NE.0) RETURN
        CALL ACTOBJ('MMODEL  ',ipmod2,1)
        IF (IERR.NE.0) RETURN
        mmode2 = ipmod2
        n2 = mmode2.kmodel(/1)
        if (n2.ne.1) then
          write(ioimp,*) 'Liaison conditionnelle mal specifiee (1)'
          call erreur(5)
          return
        endif
        imode2 = mmode2.kmodel(1)
        if (imode2.formod(1).ne.'LIAISON') THEN
          write(ioimp,*) 'Liaison conditionnelle mal specifiee (2)'
          call erreur(5)
          return
        endif
        if (klcon.eq.0) then
          nlcon = 10
          segini plicon
        endif
        klcon = klcon + 1
        if (klcon.gt.nlcon) then
          nlcon = nlcon + 10
          segadj plicon
        endif
        mlicon(klcon) = ipmod2
        tlicon(klcon) = lecon
C (fdp) option 'LIE' pour les JOI1
      ELSE IF (LECON.EQ.9) THEN
        ILIE = 1
      ELSE IF (LECON.EQ.10) THEN
        IF (LESFOR(1).NE.'MECANIQUE'.AND.LESFOR(1).NE.'POREUX') THEN
          CALL ERREUR(251)
          RETURN
        ENDIF
        CALL MODNLO(MNLOCA,NLODIM)
        IF (NLODIM.GT.NLOMAX) THEN
          CALL ERREUR(6)
          RETURN
        ENDIF
        CALL LIRMOT(MNLOCA,NLODIM,INLOC,1)
        IF (IERR.NE.0) RETURN
        CALL LIRMOT(MNLVAR,1,INLVIA,1)
        IF (IERR.NE.0) RETURN
        CALL LIROBJ('LISTMOTS',LULVIA,1,IRET)
        IF (IERR.NE.0) RETURN
      ELSE IF (LECON.GE.11.and.LECON.LE.13) THEN
        if (kbnlin.eq.0) then
          jgn = 4
          JGM = 3
          segini opnlin
        endif
        kbnlin = kbnlin + 1
        opnlin.mots(kbnlin) = mocon(lecon)
      ENDIF
      GOTO 675
 
 22   CONTINUE
      if (iptabm.gt.0.and.iptabs.eq.0.and.ipmod1.gt.0) goto 91
 
C========== ACTIVATION DU MAILLAGE POUR CONSTRUIRE MMODEL ==============
C  Recuperation des caracteristiques du MAILLAGE dans MELEME
      IF (IPGEOM .EQ. 0) THEN
        MOTERR='MAILLAGE'
        CALL ERREUR(471)
        RETURN
      ENDIF
 
      CALL ACTOBJ('MAILLAGE',IPGEOM,1)
      MELEME = IPGEOM
      NSOU  = MELEME.LISOUS(/1)
      NSOU1 = MAX(1,NSOU)
 
C=DEB==== FORMULATION HHO ==== Premieres verifications =================
      IPLHHO = 0
      IF (loHHO) THEN
C= Pour l'instant, HHO en formulation MECANIQUE !
        IF ( (NFOR.EQ.1 .AND. LESFOR(1).NE.'MECANIQUE') .OR. 
     &       (NFOR.NE.1) ) THEN
          write(ioimp,*) 'Formulation HHO --> MECANIQUE uniquement'
          CALL ERREUR(251)
          RETURN
        END IF
        IF ( .NOT. ( IFOMOD.EQ.-1 .AND. IFOUR.NE.-3) ) THEN
          write(ioimp,*) 'Formulation HHO --> 2D PLAN DEFO/CONT'
c-dbg        IF ( .NOT. ( (IFOMOD.EQ.2)                   .OR.
c-dbg     &               (IFOMOD.EQ.-1 .AND. IFOUR.NE.-3) ) ) THEN
c-dbg          write(ioimp,*) 'Formulation HHO --> 2D PLAN DEFO/CONT or 3D'
          CALL ERREUR(251)
          RETURN
        END IF
C= 
        CALL HHOPRE(CHAHHO,IPGEOM,IPLHHO,iret)
        IF (iret.NE.0) THEN
          CALL ERREUR(iret)
          RETURN
        ENDIF
      END IF
C=FIN==== FORMULATION HHO ==============================================
 
C= PARTIE 2 ============================================================
C  Initialisations et remplissage du segment MMODEL = IPMODE
C=======================================================================
c-dbg      write(ioimp,*)
c-dbg      write(ioimp,*) ' INITIALISATION MMODEL A ',NSOU1,' ZONE(S)'
 
      N1 = NSOU1
*  mmode2 sert a ranger la deuxieme partie du contact symetrique
      SEGINI,MMODEL,mmode2
      IPMODE = MMODEL
 
C  Nom du constituant par defaut si non donne en entree
      IF (kcons.EQ.0) WRITE(CONM,FMT='(I16)') IPMODE
 
C  Determination de MN3 selon la formulation :
C  Par defaut MN3 est fixe a 1.
      MN3 = 1
      IF (LESFOR(1).EQ.'NAVIER_STOKES   ') THEN
        MN3=2
C* NAVIER_STOKES + NLIN idem que MECANIQUE
        IF (LESPRO(1).EQ.'NLIN    ') MN3=12
      ENDIF
      IF (LESFOR(1).EQ.'EULER           ') MN3=2
      IF (LESFOR(1).EQ.'DARCY           ') MN3=2
      IF (LESFOR(1).EQ.'CHANGEMENT_PHASE') MN3=12
      IF (LESFOR(1).EQ.'THERMOHYDRIQUE  ') MN3=12
      IF (LESFOR(1).EQ.'METALLURGIE     ') MN3=2
      IF (LESFOR(1).EQ.'FISSURE         ') MN3=2
      IF (LESFOR(1).EQ.'MECANIQUE       ' .OR.
     &    LESFOR(1).EQ.'POREUX          ' .OR.
     &    LESFOR(1).EQ.'CHARGEMENT      ' .OR.
     &    nfor.EQ.2 ) THEN
         IF (INLOC.NE.0) THEN
           MN3=14
         ELSE
           MN3=12
         ENDIF
      ENDIF
      IF (LESFOR(1).EQ.'LIQUIDE         ') MN3=12
      IF (LESFOR(1).EQ.'LIAISON         ') MN3=12
      IF (LESFOR(1).EQ.'ELECTROSTATIQUE ') MN3=12
      IF (LESFOR(1).EQ.'DIFFUSION       ') MN3=12
      IF (LESFOR(1).EQ.'MAGNETODYNAMIQUE') MN3=12
      IF (LESFOR(1).EQ.'MELANGE         ') MN3= 7
c-dbg      write(ioimp,*) 'MN3 =',MN3,LESFOR(1),nfor
 
C- Determination de NOBMOD selon la formulation :
      NOBMOD = 0
      NOBMEC = 0
      NOBDIF = 0
      IF     (LESFOR(1).EQ.'CONTACT         ') THEN
c*        IF (IFROCA.NE.0) NOBMOD = 2
c*        IF (ifrtt .ne.0) NOBMOD = 1
        NOBMOD = 3
 
      ELSEIF     (LESFOR(1).EQ.'CONTRAINTE') THEN
        if (lactr.eq.1 .and. idim.ne.3) nobmod=3
        if (lactr.eq.1 .and. idim.eq.3) nobmod=4
        if (lactr.eq.2) nobmod=3
 
      ELSEIF (LESFOR(1).EQ.'DIFFUSION       ') THEN
        NOBDIF = NOBMOD
        NOBMOD = NOBMOD + 1
        IF (LDILOI.GT.0) NOBMOD = NOBMOD + 4
 
C* Modeles UTILISATEUR en MECANIQUE :
      ELSEIF (LMEEXT) THEN
        NOBMEC = NOBMOD
        IF (LMELOI.GT.0) NOBMOD = NOBMOD + 4
C        IF (LMEVIX     ) NOBMOD = NOBMOD + 2
        IF (LMEVIX     ) NOBMOD = NOBMOD + 1
 
      ELSEIF (LESFOR(1).eq.'METALLURGIE     ') THEN
C       On rangera les pointeurs sur les ListMots Phases,
C       Reactifs, Produits et Types de Reactions dans IVAMOD
        NOBMOD = 4
 
      ELSEIF (LESFOR(1).eq.'CHANGEMENT_PHASE') THEN
 
        IF     (LESPRO(1)(1:10).EQ.'PARFAIT   ') THEN
C       On rangera : -le LISTMOTS des inconnues primales et duales dedans
C                    -le MAILLAGE des MULTIPLICATEURS 'LX'
          NOBMOD = 2
        ELSEIF (LESPRO(1)(1:10).EQ.'SOLUBILITE') THEN
C       On rangera : -le LISTMOTS des inconnues primales et duales dedans
C                    -le 1er  MAILLAGE des MULTIPLICATEURS 'LX'
C                    -le 2eme MAILLAGE des MULTIPLICATEURS 'LX'
          NOBMOD = 3
        ELSE
          CALL ERREUR(5)
        ENDIF
 
      ELSE IF (lesfor(1).eq.'NAVIER_STOKES   ') THEN
        IF (LESPRO(1).EQ.'NLIN    ') NOBMOD = 1
 
      ELSE IF (lesfor(1).eq.'THERMIQUE       ') THEN
        IF (IRAYE.eq.1) NOBMOD = 2*icavit+isyme*idim+ifacaf*4
 
C LIAISON conditionnelle
      ELSE IF (LESFOR(1).EQ.'LIAISON         ') THEN
c*        if (klcon.gt.0) THEN
          noblia = NOBMOD
          NOBMOD = NOBMOD + klcon
c*        endif
      ENDIF
 
C=DEB==== FORMULATION HHO ==== Donnees supplementaires =================
      IF (loHHO) THEN
        nobHHO = NOBMOD
        NOBMOD = NOBMOD + MTYHHO
      END IF
C=FIN==== FORMULATION HHO ==============================================
 
      IF (jlmot1.gt.0) THEN
        if (nobmod.ne.0) then
          write(ioimp,*) 'NOBMOD !=0 et jlmot1,jlmmot2 > 0'
          call erreur(5)
          return
        endif
        NOBMOD = 2
      ENDIF
 
      NOBMOD0 = NOBMOD
      NOBMEC0 = NOBMEC
      NOBDIF0 = NOBDIF
      nobHHO0 = nobHHO
      noblia0 = noblia
 
C***********************************************************************
C  Boucle (10) sur les maillages elementaires de IPGEOM
C***********************************************************************
      MELEME = IPGEOM
      IPT1 = MELEME
 
      DO 10 IM = 1, NSOU1
 
        IF (NSOU.NE.0) IPT1 = MELEME.LISOUS(IM)
        ITYP1 = IPT1.ITYPEL
        NBNN  = IPT1.NUM(/1)
        NBEL  = IPT1.NUM(/2)
 
C +--------------------------------------------------------------------+
C | Creation du modele elementaire IMODEL                              |
C +--------------------------------------------------------------------+
        NOBMOD = NOBMOD0
        NOBMEC = NOBMEC0
        NOBDIF = NOBDIF0
        nobHHO = nobHHO0
        noblia = noblia0
 
        SEGINI,IMODEL
        mmodel.KMODEL(IM) = IMODEL
 
C +--------------------------------------------------------------------+
C | Remplissage du IMODEL                                              |
C +--------------------------------------------------------------------+
        imodel.IMAMOD = IPT1
 
        imodel.CONMOD(1:16)  = CONM
        imodel.conmod(17:24) = PHAM
 
        DO i = 1, NFOR
          imodel.FORMOD(i) = LESFOR(i)
        ENDDO
        IF (NMAT.NE.0) THEN
          DO i = 1, NMAT
            imodel.MATMOD(i) = LESPRO(i)
          ENDDO
        ENDIF
 
C  Informations liees au MATERIAU/COMPORTEMENT
        CMATE = '        '
        IMATE = 0
        INATU = 0
        CALL NOMATE(FORMOD,NFOR,MATMOD,NMAT,CMATE,IMATE,INATU)
c-dbg      write(6,*)'NOMATE : ',cmate,imate,inatu
c-dbg      write(6,*)'  MATMOD =',(matmod(i),i=1,nmat)
        IF (IERR.NE.0) then
          write(ioimp,*) ' Probleme apres NOMATE'
          CALL ERREUR(251)
          GOTO 990
        ENDIF
C* Petit cas particulier en cas de modele VISCO_EXTERNE :
C* On recupere IVIEX stocke dans INATU (cf. NOMATE)
        IF (INATU .LE. -2) THEN
          IF (.NOT. LMEVIX) THEN
      write(ioimp,*)'  VISCO_EXTERNE : INATU & LMEVIX incompatibles'
            call erreur(5)
            return
          ENDIF
          IVIEX = -2 - INATU
          INATU = -2
        ENDIF
        imodel.CMATEE = CMATE
        imodel.IMATEE = IMATE
        imodel.INATUU = INATU
 
c        imodel.ideriv = jderiv
cbp,2020-12-10 : abandon de MEPSIL (CCOPTIO) et IDERIV (MMODEL)
        imodel.IDERIV = 0
 
        ipmmel = 0
        IF (LESFOR(1).EQ.'MELANGE         ') THEN
          IF (CMATE.EQ.'PARALLEL' .OR. CMATE.EQ.'SERIE') THEN
            ipmmel = LESMOD(1)
            if (ipmmel.le.0) then
              call erreur(21)
              GOTO 990
            endif
          ENDIF
        ENDIF
 
C +--------------------------------------------------------------------+
C | Remplissage des couples TYMODE/IVAMOD                              |
C +--------------------------------------------------------------------+
        IF (LESFOR(1).EQ.'CONTACT         ') THEN
          IF (IFROCA.EQ.0) THEN
*  lecture du mot-cle
            call lirmot(MCTCT,4,iret,0)
            ictct=iret
            if(iret.eq.0) ictct=1
            ictr=2
            if(ictct.ne.2) ictr=1
*  cas mortar : uniquement disponible en 2D
            if (iret.eq.4) then
              if (idim.ne.2) then
                INTERR(1) = IDIM
                CALL ERREUR(1104)
                GOTO 990
              endif
              ictr=4
            endif
*  lecture du deuxieme maillage
            call lirobj('MAILLAGE',ipgeo2,1,iretou)
            if (ierr.ne.0) return
            ipgeox=ipgeo2
            call mocon1(ipgeox,lecont,ictr)
            if(ierr.ne.0) return
            tymode(1)='MAILLAGE'
            ivamod(1)=ipgeom
            tymode(2)='MAILLAGE'
            ivamod(2)=ipgeo2
            tymode(3)='ENTIER'
            ivamod(3)=ictr
            imamod=ipgeox
            if(ictct.eq.2) then
              ip2=ipgeom
              call mocon1(ip2,lecont,ictr)
              if(ierr.ne.0) return
              ip1=imamod
              call fuse(ip1,ip2,iret,.false.)
              if(ierr.ne.0) return
              imamod=iret
            endif
            if(ictct.eq.3) then
              segini,imode1
              mmode2.kmodel(im)=imode1
              ipgeox=ipgeom
              call mocon1(ipgeox,lecont,ictr)
              if(ierr.ne.0) return
              imode1.tymode(1)='MAILLAGE'
              imode1.ivamod(1)=ipgeo2
              imode1.tymode(2)='MAILLAGE'
              imode1.ivamod(2)=ipgeom
              imode1.tymode(3)='ENTIER'
              imode1.ivamod(3)=1
              imode1.imamod=ipgeox
            endif
          ENDIF
 
          IF (IFROCA.EQ.1) THEN
*  deuxieme maillage deja lu
            ipgeo2 = ipgeom
            ictr=0
            call mocon1(ipgeo2,lecont,ictr)
            if(ierr.ne.0) return
            imamod=ipgeo2
            ipt3=ipgeo2
            segact ipt3
            ityp1=ipt3.itypel
            TYMODE(1)='MAILLAGE'
            IVAMOD(1)=IPGEOM
            TYMODE(2)='MAILLAGE'
            IVAMOD(2)=IBETON
          ENDIF
**          IF (ifrtt.eq.1) then
**            ivamod(1)=ipgeo2
**            tymode(1)='MAILLAGE'
**            if (ipgeo2.eq.0) then
**             call erreur(641)
**             return
**            endif
**          ENDIF
 
        ELSEIF (lesfor(1).eq.'CONTRAINTE      ') then
          if (lactr.eq.1.or.lactr.eq.2)  then
            call mocon2(ipgeom,ipt7)
          endif
          if (lactr.eq.3) call mocon3(ipgeom,ipt7)
          tymode(1)='ENTIER'
          ivamod(1)=lactr
          tymode(2)='MAILLAGE'
          ivamod(2)=ipgeom 
          imamod=ipt7 
*  cas rotation idim-1 pts
          if (lactr.eq.1) then
            call meslir(0)
            call lirobj('POINT',ip1,1,iok)
            if (idim.eq.3) call lirobj('POINT',ip2,1,iok)
            if (ierr.ne.0) return
            tymode(3)='POINT'
            ivamod(3)=ip1    
            if (idim.eq.3) then
              tymode(4)='POINT'
              ivamod(4)=ip2    
            endif
*  cas deplacement  1 pt
          elseif (lactr.eq.2) then
            call lirobj('POINT',ip1,1,iok)
            if (ierr.ne.0) return
            tymode(3)='POINT'
            ivamod(3)=ip1    
          endif
 
        ELSEIF (lesfor(1).eq.'NAVIER_STOKES   ') THEN
          IF (nobmod.gt.0) THEN
            tymode(1) = 'LISTMOTS'
            ivamod(1) = opnlin
          ENDIF
 
        ELSEIF (lesfor(1).eq.'METALLURGIE     ') then
C         lucvar : les noms des phases
          IVAMOD(1) = lucvar
          TYMODE(1) = 'LISTMOTS'
C           reacti : les noms des reactifs
          IVAMOD(2) = ireact
          TYMODE(2) = 'LISTMOTS'
C           produi : les noms des produits
          IVAMOD(3) = iprodu
          TYMODE(3) = 'LISTMOTS'
C           lucmat : les noms des types de reactions
          IVAMOD(4) = lucmat
          TYMODE(4) = 'LISTMOTS'
 
        ELSEIF (lesfor(1).eq.'CHANGEMENT_PHASE') then
C           ipridu : les noms des variables primales et duales
          IVAMOD(1) = ipridu
          TYMODE(1) ='LISTMOTS'
          CALL IMPP1(IPT1,ipgeo2,ipgeo3,LESPRO(1))
C           ipgeo2 & ipgeo3 : MAILLAGE support des Multiplicateurs de Lagrange ('MULT')
          IF     (LESPRO(1)(1:10).EQ.'PARFAIT   ') THEN
            IVAMOD(2) = ipgeo2
            TYMODE(2) ='MAILLAGE'
          ELSEIF (LESPRO(1)(1:10).EQ.'SOLUBILITE') THEN
            IVAMOD(2) = ipgeo2
            TYMODE(2) ='MAILLAGE'
            IVAMOD(3) = ipgeo3
            TYMODE(3) ='MAILLAGE'
          ELSE
            CALL ERREUR(5)
          ENDIF
 
        ELSEIF (LESFOR(1).EQ.'THERMIQUE')THEN
          if(IRAYE.ne.0) then
            if(icavit.ne.0) then
              tymode(1)='ENTIER'
              ivamod(1)=nbga
              tymode(2)='ENTIER'
              ivamod(2)=nbdang
              if(isyme.eq.1) then
                tymode(3)='POINT'
                tymode(4)='POINT'
                if(idim.eq.3)tymode(5)='POINT'
                ivamod(3)=ipp1
                ivamod(4)=ipp2
                if(idim.eq.3)ivamod(5)=ipp3
              endif
            endif
            if(ifacaf.ne.0) then
              tymode(1)='MAILLAGE'
              tymode(2)='MAILLAGE'
              tymode(3)='MAILLAGE'
              tymode(4)='MMODEL'
              ivamod(1)= ipfac1
              ivamod(2)= ipfac2
              ivamod(3)= ipfac3
              ivamod(4)= imoco
            endif
          endif
 
        ELSEIF (LESFOR(1).EQ.'DIFFUSION') THEN
          JGN = LOCOMP
          JGM = 2
          SEGINI,MLMOT1
          TYMODE(NOBDIF+1)='LISTMOTS'
          IVAMOD(NOBDIF+1)=iplrdi
          NOBDIF = NOBDIF+1
          IF (LDILOI.GT.0) THEN
C           Indicateur 'LDIEXT' pour retrouver ses petits
            CALL POSCHA('LDIEXT  ',I_POS)
            TYMODE(NOBDIF+1)='MOT     '
            IVAMOD(NOBDIF+1)= I_POS
 
C           Pointeur vers la loi (donne par PTRLOI)
            TYMODE(NOBDIF+2)='ENTIER  '
            IVAMOD(NOBDIF+2)= LDIPTR
 
C           LMELIB : Nom de la bibliotheque (sans chemin et extension)
            CALL POSCHA(LDILIB(1:LDILGB),I_POS)
            TYMODE(NOBDIF+3)='MOT     '
            IVAMOD(NOBDIF+3)= I_POS
 
C           LMEFCT : Nom de la fonction (dans la bibliotheque)
            CALL POSCHA(LDIFCT(1:LDILGT),I_POS)
            TYMODE(NOBDIF+4)='MOT     '
            IVAMOD(NOBDIF+4)= I_POS
            NOBDIF = NOBDIF + 4
          ENDIF
 
        ELSEIF (LMEEXT) THEN
C           Modeles utilisateur en MECANIQUE :
          IF (LMELOI.GT.0) THEN
C           Indicateur 'LMEEXT' pour retrouver ses petits
            CALL POSCHA('LMEEXT  ',I_POS)
            TYMODE(NOBMEC+1)='MOT     '
            IVAMOD(NOBMEC+1)= I_POS
 
C           Pointeur vers la loi (donne par PTRLOI)
            TYMODE(NOBMEC+2)='ENTIER  '
            IVAMOD(NOBMEC+2)= LMEPTR
 
C           LMELIB : Nom de la bibliotheque (sans chemin et extension)
            CALL POSCHA(LMELIB(1:LMELGB),I_POS)
            TYMODE(NOBMEC+3)='MOT     '
            IVAMOD(NOBMEC+3)= I_POS
 
C           LMEFCT : Nom de la fonction (dans la bibliotheque)
            CALL POSCHA(LMEFCT(1:LMELGT),I_POS)
            TYMODE(NOBMEC+4)='MOT     '
            IVAMOD(NOBMEC+4)= I_POS
            NOBMEC = NOBMEC + 4
          ENDIF
 
          IF (LMEVIX) THEN
c*          IF (INATU .EQ. -2) THEN
            NOBMEC = NOBMEC + 1
            imodel.TYMODE(NOBMEC) = 'IVIEX   '
            imodel.IVAMOD(NOBMEC) = IVIEX
C            NOBMEC = NOBMEC + 1
C            imodel.TYMODE(NOBMEC) = 'ENTIER  '
C            imodel.IVAMOD(NOBMEC) = 0
          ENDIF
 
        ELSE IF (LESFOR(1).EQ.'LIAISON         ') THEN
          if (klcon.gt.0) THEN
            do i = 1, klcon
              if (tlicon(i).eq.6) TYMODE(noblia+i) = 'CONDINFE'
              if (tlicon(i).eq.7) TYMODE(noblia+i) = 'CONDSUPE'
              IVAMOD(noblia+i) = mlicon(i)
            enddo
          ENDIF
 
        ENDIF
        if (jlmot1.gt.0) then
          ivamod(1) = jlmot1
          ivamod(2) = jlmot2
          tymode(1) = 'LISTMOTS'
          tymode(2) = 'LISTMOTS'
          lucvar = jlmot1
          lucmat = jlmot2
        endif
 
C=DEB==== FORMULATION HHO ==== Remplissage de donnees ==================
        IF (loHHO) THEN
          modHHO = imodel
          CALL HHOPRM(chaHHO,modHHO,nobHHO,iplHHO,iret)
          IF (iret.NE.0) THEN
            CALL ERREUR(iret)
            GOTO 990
          END IF
        END IF
C=FIN==== FORMULATION HHO ==============================================
 
c*        IF (FORMOD(1).EQ.'MELANGE         ') THEN
        IF (ipmmel.GT.0) THEN
          mmode1 = ipmmel
          n1mel = mmode1.kmodel(/1)
          NOBMOI = ivamod(/1)
          NOBMOD = NOBMOI + n1mel
          SEGADJ,IMODEL
          kbmod = 0
          DO i = 1, n1mel
            imode1 = mmode1.kmodel(i)
            if (imode1.imamod.eq.imamod) then
              if (kbmod.eq.0) then
                imode2 = imode1
              else
                if (imode1.formod(1).ne.imode2.formod(1).or.
     &              imode1.imatee.ne.imode2.imatee) goto 117
              endif
              kbmod = kbmod + 1
              tymode(NOBMOI+kbmod) = 'IMODEL'
              ivamod(NOBMOI+kbmod) = imode1
            endif
 117        continue
          ENDDO
          if (kbmod.eq.0) then
            call erreur(21)
            return
          endif
          if (kbmod.ne.n1mel) then
            NOBMOD = NOBMOI + kbmod
            SEGADJ,imodel
          endif
        ENDIF
c*        ENDIF
c-dbg      write(ioimp,*) 'Fin remplissage IVAMOD/TYMODE'
 
C +--------------------------------------------------------------------+
C | Determination de la valeur de NEFMOD pour IMODEL                   |
C +--------------------------------------------------------------------+
C  Affectation du type d'ELEMENTS FINIS si donnes par utilisateur
 
C     Cas des SURE (relation de conformite) :
C       NEPAPA = si EF specifique demande -> on utilise ses inconnues
        NEPAPA = 0
        IF (ITYP1.EQ.48) THEN
          imodel.NEFMOD = 259
          IF (ITEF.GT.0) THEN
            DO i=1,ITEF
              CALL PLACE(NOMTP,LNOMTP,MELE,LESTEF(i))
              IF (MELE.NE.0) NEPAPA = MELE
            ENDDO
          ENDIF
          IF (NEPAPA.EQ.0) THEN
c         2D -> on choisit les inconnues du QUA4 pour toute formulation
            IF (IDIM.EQ.2) THEN
              NEPAPA=8
c         3D -> on choisit les inconnues du CUB8 pour toute formulation
            ELSEIF (IDIM.EQ.3) THEN
              NEPAPA=14
            ELSE
              CALL ERREUR(610)
              RETURN
            ENDIF
          ENDIF
          MELE = NEPAPA
          GOTO 101
        ENDIF
 
C=DEB==== FORMULATION HHO ==== NEFMOD = HHO_NUM_ELEMENT pour tous les elements =====
        IF (loHHO) THEN
          imodel.NEFMOD = HHO_NUM_ELEMENT
          MELE = imodel.NEFMOD
          GOTO 101
        END IF
C=FIN==== FORMULATION HHO ==============================================
 
        NEFMOD = 0
        IF (ITEF.NE.0) THEN
C  Cas de la FORMULATION 'NAVIER_STOKES'
          IF (LESFOR(1).EQ.'NAVIER_STOKES') THEN
            IF (LESTEF(1).EQ.'LINE')THEN
              IF (ITYP1.EQ. 3) NEFMOD=129
              IF (ITYP1.EQ. 7) NEFMOD=130
              IF (ITYP1.EQ.11) NEFMOD=131
              IF (ITYP1.EQ.33) NEFMOD=132
              IF (ITYP1.EQ.34) NEFMOD=133
              IF (ITYP1.EQ.35) NEFMOD=134
              IF (ITYP1.EQ.36) NEFMOD=135
            ELSE IF(LESTEF(1).EQ.'MACR')THEN
              IF (ITYP1.EQ. 3) NEFMOD=136
              IF (ITYP1.EQ. 7) NEFMOD=137
              IF (ITYP1.EQ.11) NEFMOD=138
              IF (ITYP1.EQ.33) NEFMOD=139
              IF (ITYP1.EQ.34) NEFMOD=140
              IF (ITYP1.EQ.35) NEFMOD=141
              IF (ITYP1.EQ.36) NEFMOD=142
C  Il nous manque la pyramide
            ELSE IF (LESTEF(1).EQ.'QUAF') THEN
              IF (ITYP1.EQ. 3) NEFMOD=143
              IF (ITYP1.EQ. 7) NEFMOD=144
              IF (ITYP1.EQ.11) NEFMOD=145
              IF (ITYP1.EQ.33) NEFMOD=146
              IF (ITYP1.EQ.34) NEFMOD=147
              IF (ITYP1.EQ.35) NEFMOD=148
              IF (ITYP1.EQ.36) NEFMOD=149
C Il nous manque la pyramide
            ELSE IF (LESTEF(1).EQ.'LINB') THEN
              IF (ITYP1.EQ. 3) NEFMOD=158
              IF (ITYP1.EQ. 7) NEFMOD=159
              IF (ITYP1.EQ.11) NEFMOD=160
              IF (ITYP1.EQ.33) NEFMOD=161
              IF (ITYP1.EQ.34) NEFMOD=162
C             IF (ITYP1.EQ.35) NEFMOD=163
C             IF (ITYP1.EQ.36) NEFMOD=164
C Il nous manque la pyramide et le tetrahedre
            ELSE
              DO i=1,ITEF
                CALL PLACE(NOMTP,LNOMTP,MELE,LESTEF(i))
                IF (MELE.EQ.0) GOTO 99
                MEGE=NUMGEO(MELE)
                IF (MEGE.EQ.0) GOTO 99
                IF (MEGE.EQ.ITYP1) GOTO 610
              ENDDO
              GOTO 99
 610          NEFMOD=MELE
            ENDIF
C  Cas de la FORMULATION 'EULER'
          ELSE IF (LESFOR(1).EQ.'EULER') THEN
            IF (ITYP1.EQ. 2) NEFMOD=ITYP1
            IF (ITYP1.EQ. 4) NEFMOD=ITYP1
            IF (ITYP1.EQ. 8) NEFMOD=ITYP1
            IF (ITYP1.EQ.14) NEFMOD=ITYP1
            IF (ITYP1.EQ.16) NEFMOD=ITYP1
            IF (ITYP1.EQ.23) NEFMOD=ITYP1
            IF (ITYP1.EQ.25) NEFMOD=ITYP1
C  Cas des autres FORMULATIONs
          ELSE
            DO i=1,ITEF
              if(lestef(i)(1:4).eq.'BBAR') lobbar = .true.
              if (lobbar) CALL MODE20(ITYP1,LESTEF(I))
              CALL PLACE(NOMTP,LNOMTP,MELE,LESTEF(i))
              IF (MELE.EQ.0) GOTO 99
              MEGE=NUMGEO(MELE)
              IF (MEGE.EQ.0) GOTO 99
              IF (MEGE.EQ.ITYP1) GOTO 6
c kich cas du POI1
              if (ityp1.eq.1) goto 6
            ENDDO
            GOTO 99
C  Cas particulier pour les elements polygonaux
 6          IF (ITYP1.EQ.32) THEN
              MELE=MELE+NBNN-3
              IF (NBNN.GT.14) GOTO 99
            ENDIF
            NEFMOD=MELE
          ENDIF
C  Affectation des elements finis de maniere automatique
        ELSE
C  Cas des milieux POREUX
          IF (LESFOR(1).EQ.'POREUX') THEN
            IF (ITYP1.EQ. 6) NEFMOD=79
            IF (ITYP1.EQ.10) NEFMOD=80
            IF (ITYP1.EQ.15) NEFMOD=81
            IF (ITYP1.EQ.24) NEFMOD=82
            IF (ITYP1.EQ.17) NEFMOD=83
            IF (ITYP1.EQ.29) NEFMOD=108
            IF (ITYP1.EQ.30) NEFMOD=109
            IF (ITYP1.EQ.31) NEFMOD=110
C  Cas des elements de frottement (formulation FROTTEMENT)
          ELSE IF (LESFOR(1).EQ.'CONTACT') THEN
            NEFMOD=22                                    
            if(ifrtt.eq.1) then
              IF (ITYP1.EQ.22.AND.IDIM.EQ.2) NEFMOD=107
              IF (ITYP1.EQ.22.AND.IDIM.EQ.3) NEFMOD=165
            elseif(ifroca.ne.0) then
              IF (ITYP1.EQ.22.AND.IDIM.EQ.2)  NEFMOD=261
              IF (ITYP1.EQ.22.AND.IDIM.EQ.3)  NEFMOD=262
            endif
C  Cas des elements de contrainte (formulation CONTRAINTE)
          ELSE IF (LESFOR(1).EQ.'CONTRAINTE') THEN
             NEFMOD=22
C  Cas des elements hybrides (imposes en DARCY)
          ELSE IF (LESFOR(1).EQ.'DARCY') THEN
            IF (ITYP1.EQ. 3) NEFMOD=143
C           IF (ITYP1.EQ. 4) NEFMOD=99
C           IF (ITYP1.EQ. 8) NEFMOD=100
C           IF (ITYP1.EQ.23) NEFMOD=101
C           IF (ITYP1.EQ.16) NEFMOD=102
C           IF (ITYP1.EQ.14) NEFMOD=103
            IF (ITYP1.EQ. 7) NEFMOD=99
            IF (ITYP1.EQ.11) NEFMOD=100
            IF (ITYP1.EQ.35) NEFMOD=101
            IF (ITYP1.EQ.34) NEFMOD=102
            IF (ITYP1.EQ.33) NEFMOD=103
C  Cas de la formulation MAGNETODYNAMIQUE
          ELSE IF (LESFOR(1).EQ.'MAGNETODYNAMIQUE') THEN
            IF (ITYP1.EQ. 4) NEFMOD=128
C  Cas de la formulation 'NAVIER_STOKES'
          ELSE IF (LESFOR(1).EQ.'NAVIER_STOKES') THEN
            IF (ILNAVI.EQ.0) THEN
              CALL MESLIR(-341)
              GOTO 990
            ELSEIF (ILNAVI.EQ.1) THEN
C         LICE
              IF (ITYP1.EQ. 3) NEFMOD=195
              IF (ITYP1.EQ. 7) NEFMOD=196
              IF (ITYP1.EQ.11) NEFMOD=197
              IF (ITYP1.EQ.33) NEFMOD=198
              IF (ITYP1.EQ.34) NEFMOD=199
              IF (ITYP1.EQ.35) NEFMOD=200
              IF (ITYP1.EQ.36) NEFMOD=201
            ELSEIF (ILNAVI.EQ.2) THEN
C         LIMS
              IF (ITYP1.EQ. 3) NEFMOD=202
              IF (ITYP1.EQ. 7) NEFMOD=203
              IF (ITYP1.EQ.11) NEFMOD=204
              IF (ITYP1.EQ.33) NEFMOD=205
              IF (ITYP1.EQ.34) NEFMOD=206
              IF (ITYP1.EQ.35) NEFMOD=207
              IF (ITYP1.EQ.36) NEFMOD=208
            ELSEIF (ILNAVI.EQ.3) THEN
C         LBMS
              IF (ITYP1.EQ. 3) NEFMOD=209
              IF (ITYP1.EQ. 7) NEFMOD=210
              IF (ITYP1.EQ.11) NEFMOD=211
              IF (ITYP1.EQ.33) NEFMOD=212
              IF (ITYP1.EQ.34) NEFMOD=213
              IF (ITYP1.EQ.35) NEFMOD=214
              IF (ITYP1.EQ.36) NEFMOD=215
            ELSEIF (ILNAVI.EQ.4) THEN
C         MCCE
              IF (ITYP1.EQ. 3) NEFMOD=216
              IF (ITYP1.EQ. 7) NEFMOD=217
              IF (ITYP1.EQ.11) NEFMOD=218
              IF (ITYP1.EQ.33) NEFMOD=219
              IF (ITYP1.EQ.34) NEFMOD=220
              IF (ITYP1.EQ.35) NEFMOD=221
              IF (ITYP1.EQ.36) NEFMOD=222
            ELSEIF (ILNAVI.EQ.5) THEN
C         MCP1
              IF (ITYP1.EQ. 3) NEFMOD=223
              IF (ITYP1.EQ. 7) NEFMOD=224
              IF (ITYP1.EQ.11) NEFMOD=225
              IF (ITYP1.EQ.33) NEFMOD=226
              IF (ITYP1.EQ.34) NEFMOD=227
              IF (ITYP1.EQ.35) NEFMOD=228
              IF (ITYP1.EQ.36) NEFMOD=229
            ELSEIF (ILNAVI.EQ.6) THEN
C         MCMS
              IF (ITYP1.EQ. 3) NEFMOD=230
              IF (ITYP1.EQ. 7) NEFMOD=231
              IF (ITYP1.EQ.11) NEFMOD=232
              IF (ITYP1.EQ.33) NEFMOD=233
              IF (ITYP1.EQ.34) NEFMOD=234
              IF (ITYP1.EQ.35) NEFMOD=235
              IF (ITYP1.EQ.36) NEFMOD=236
            ELSEIF (ILNAVI.EQ.7) THEN
C         QFCE
              IF (ITYP1.EQ. 3) NEFMOD=237
              IF (ITYP1.EQ. 7) NEFMOD=238
              IF (ITYP1.EQ.11) NEFMOD=239
              IF (ITYP1.EQ.33) NEFMOD=240
              IF (ITYP1.EQ.34) NEFMOD=241
              IF (ITYP1.EQ.35) NEFMOD=242
              IF (ITYP1.EQ.36) NEFMOD=243
            ELSEIF (ILNAVI.EQ.8) THEN
C         QFP1
              IF (ITYP1.EQ. 3) NEFMOD=244
              IF (ITYP1.EQ. 7) NEFMOD=245
              IF (ITYP1.EQ.11) NEFMOD=246
              IF (ITYP1.EQ.33) NEFMOD=247
              IF (ITYP1.EQ.34) NEFMOD=248
              IF (ITYP1.EQ.35) NEFMOD=249
              IF (ITYP1.EQ.36) NEFMOD=250
            ELSEIF (ILNAVI.EQ.9) THEN
C         QFMS
              IF (ITYP1.EQ. 3) NEFMOD=251
              IF (ITYP1.EQ. 7) NEFMOD=252
              IF (ITYP1.EQ.11) NEFMOD=253
              IF (ITYP1.EQ.33) NEFMOD=254
              IF (ITYP1.EQ.34) NEFMOD=255
              IF (ITYP1.EQ.35) NEFMOD=256
              IF (ITYP1.EQ.36) NEFMOD=257
            ENDIF
C  Cas de la formulation 'EULER'
          ELSE IF (LESFOR(1).EQ.'EULER') THEN
            IF (ITYP1.EQ. 2) NEFMOD=ITYP1
            IF (ITYP1.EQ. 4) NEFMOD=ITYP1
            IF (ITYP1.EQ. 8) NEFMOD=ITYP1
            IF (ITYP1.EQ.14) NEFMOD=ITYP1
            IF (ITYP1.EQ.16) NEFMOD=ITYP1
            IF (ITYP1.EQ.23) NEFMOD=ITYP1
            IF (ITYP1.EQ.25) NEFMOD=ITYP1
C  Cas des autres formulations
          ELSE
            NEFMOD=ITYP1
c kich cas du POI1
            if (ityp1.eq.1) nefmod = 45
C  Cas particuliers des elements polygonaux
            IF (NEFMOD.EQ.32) NEFMOD=111+NBNN-3
c gounand cas des 'CU27','PR21','TE15','PY19'
            if (NEFMOD.GE.33.AND.NEFMOD.LE.36) then
               nefmod = nefmod-33 +275
            endif
C  Cas particuliers des elements finis pour IDIM=1
            IF (IDIM.EQ.1) THEN
              NEFMOD = 0
              IF (LESFOR(1).EQ.'THERMIQUE') THEN
                IF (ICONV.NE.0 .OR. IRAYE.NE.0) THEN
                  IF (ITYP1.EQ.1) NEFMOD=45
                  IF (ITYP1.EQ.2) NEFMOD=ITYP1
                ELSE
                  IF (ITYP1.EQ.2) NEFMOD=191
                  IF (ITYP1.EQ.3) NEFMOD=192
                ENDIF
              ELSE IF (LESFOR(1).EQ.'MECANIQUE') THEN
                IF (ITYP1.EQ.2) NEFMOD=193
                IF (ITYP1.EQ.3) NEFMOD=194
              ELSE IF (LESFOR(1).EQ.'FISSURE') THEN
                IF (ITYP1.EQ.2) NEFMOD=ITYP1
              ELSE IF (LESFOR(1).EQ.'ELECTROSTATIQUE') THEN
                IF (ITYP1.EQ.2) NEFMOD=193
                IF (ITYP1.EQ.3) NEFMOD=194
              ELSE IF (LESFOR(1).EQ.'DIFFUSION') THEN
* En attendant le retour a la normale pour la diffusion, on ajoute une
* enieme rustine en mettant les memes elements qu'en thermique.
**                IF (ITYP1.EQ.2) NEFMOD=193
**                IF (ITYP1.EQ.3) NEFMOD=194
                IF (ITYP1.EQ.2) NEFMOD=191
                IF (ITYP1.EQ.3) NEFMOD=192
              ENDIF
            ENDIF
          ENDIF
        ENDIF
        IF (NEFMOD.EQ.0) GOTO 99
        MELE = NEFMOD
 101    CONTINUE
C +--------------------------------------------------------------------+
C | Fin de la valeur de NEFMOD pour IMODEL                             |
C +--------------------------------------------------------------------+
      MFR  = NUMMFR(NEFMOD)
      mfr2 = NUMMFR(MELE)
c-dbg      write(6,*)'  ITYP1  =',ityp1,nefmod,mele,MFR,mfr2
 
C +--------------------------------------------------------------------+
C  Quelques tests supplementaires en attendant mieux
        IF (LESFOR(1).EQ.'THERMIQUE') THEN
C          nnz = MATMOD(/2)
          iplaz = 0
          call place(MATMOD,MATMOD(/2),iplaz,'PHASE')
          IF (iplaz.ne.0 ) THEN
c    test que les elements sont lineaires
c* on a aussi : ipt4 = ipt1 = imodel.imamod
            ipt4 = imodel.imamod
            itt = ipt4.itypel
            if (kdegre(itt) .gt. 2) then
              call erreur(982)
              goto 990
            endif
          ENDIF
        ENDIF
 
        IF (LESFOR(1).EQ.'MECANIQUE') THEN
C  Cas du materiau unidirectionnel
          IF (IMATE.EQ.4) THEN
C  Cas des cerces : sans interet !
            IF (MFR.EQ.27) THEN
              CALL ERREUR(251)
              GOTO 990
            ENDIF
C  Cas de la plasticite
            IF (INATU.NE.0) THEN
C  OK si massif bidim ou si coque tridim dans le cas acier_uni
              IF (INATU.EQ.40)THEN
                IF ((MFR.NE.1.OR.IFOUR.GT.0).AND.
     &             ((MFR.NE.3.AND.MFR.NE.9).OR.IFOUR.NE.2)) THEN
                  CALL ERREUR(251)
                  GOTO 990
                ENDIF
C  Dans les autres cas, on n'autorise pour le moment que COQ2 et massif
              ELSE IF (MELE.NE.44.AND.MFR.NE.1) THEN
                CALL ERREUR(251)
                GOTO 990
              ENDIF
            ENDIF
          ENDIF
C
C   Cas du materiau 'ZONE_COHESIVE'
          IF (IMATE.EQ.12) THEN
            IF (MFR.NE.77) THEN
              CALL ERREUR(251)
              GOTO 990
            ENDIF
          ENDIF
 
C  Cas du modele section : on n'autorise pour le moment que TIMO
          IF (CMATE.EQ.'SECTION'.AND.MELE.NE.84) THEN
            CALL ERREUR(251)
            GOTO 990
          ENDIF
        ENDIF
C  Le modele de GURSON n'est possible qu'en 3D, axisymetrique ou
C  deformations planes
        IF (INATU.EQ.38) THEN
          IF ( IFOUR.NE.0 .AND. IFOUR.NE.2 .AND. IFOUR.NE.-1 ) THEN
            MOTERR(1:8)='GURSON'
            MOTERR(9:16)='MECANIQU'
            INTERR(1) = IFOUR
            CALL ERREUR (81)
            GOTO 990
          ENDIF
        ENDIF
 
C  Le modele ISS_GRANGE n'est utilisable qu'en 3D
        IF (INATU.EQ.151 .AND. IFOUR.NE.2) THEN
          INTERR(1) = IFOUR
          CALL ERREUR (709)
          GOTO 990
        ENDIF
C  Le modele RUP_THER n'est utilisable qu'en 3D
        IF (INATU.EQ.152 .AND. IFOUR.NE.2) THEN
          INTERR(1) = IFOUR
          CALL ERREUR (709)
          GOTO 990
        ENDIF
C  Le modele COULOMB n'est utilisable qu'en 3D avec les éléments JOI1
        IF (INATU.EQ.34 .AND. IFOUR.NE.2 .AND. MFR.EQ.75) THEN
          INTERR(1) = IFOUR
          CALL ERREUR (709)
          GOTO 990
        ENDIF
C.. Restrictions en formulation 'MECANIQUE' avec une loi de
C   comportement non lineaire externe
C   Rappel : LMEEXT exprime la condition (NFOR.EQ.1) ET
C     (LESFOR(1).EQ.'MECANIQUE') ET (loi non lineaire externe)
        IF ( LMEEXT ) THEN
C         En formulation 'MECANIQUE', les lois non lineaires externes
C         n'autorisent qu'une seule composante de temperature
C         => incompatibilite avec des modeles de coques n'ayant pas
C         de points d'integration dans l'epaisseur (trois composantes
C         dans ce cas, 'TINF', 'T   ' et 'TSUP')
C         Le test ci-dessous est coherent avec celui de IDTEMP.
          IF ( (MFR.EQ.3.OR.MFR.EQ.5.OR.MFR.EQ.9).AND.
     &         (NGINT.EQ.0) ) THEN
            CALL ERREUR(951)
            GOTO 990
          ENDIF
C         Les lois de la famille 'VISCO_EXTERNE' ne s'appliquent pour
C         l'instant qu'aux elements massifs, avec option de calcul 3D
C         Et restriction pour l'instant a 'VISCO_EXTERNE' 'GENERAL'
          IF ( LMEVIX ) THEN
            IF ((MFR.NE.1.AND.MFR.NE.31) .OR. IFOUR.NE.2) THEN
              KERRE = 950
            ELSE IF ( IVIEX.NE.1 ) THEN
              KERRE = 958
            ELSE
              KERRE = 0
            ENDIF
            IF (KERRE.NE.0) THEN
              CALL ERREUR(KERRE)
              GOTO 990
            ENDIF
          ENDIF
        ENDIF
 
C  Formulation 'THERMIQUE' 'CONVECTION'
C  Adequation EF de type COQue et mot 'INFERIEURE' / 'SUPERIEURE'
        IF (ICONV.EQ.1) THEN
          CALL PLACE(LESPRO,NMAT,ISUP,'SUPERIEURE')
          CALL PLACE(LESPRO,NMAT,IINF,'INFERIEURE')
          ITOT = ISUP+IINF
          LOGRE = MELE.EQ.27 .OR. MELE.EQ.41 .OR. MELE.EQ.44 .OR.
     &            MELE.EQ.49 .OR. MELE.EQ.56
          IF ( ITOT.NE.0.AND.(.NOT. LOGRE) ) THEN
            CALL ERREUR(16)
            GOTO 990
          ENDIF
          IF ( ITOT.EQ.0 .AND. LOGRE) THEN
            CALL ERREUR(513)
            GOTO 990
          ENDIF
        ENDIF
 
C Formulation 'DIFFUSION' :
        IF (LESFOR(1) .EQ. 'DIFFUSION') THEN
C -     Verification sur les types de FORMULATION et/ou d'elements
          IF (IFOUR.EQ.2 .AND. MELE.GE.4 .AND. MELE.LT.11) THEN
            CALL ERREUR(16)
            GOTO 99
          ENDIF
          IF (MFR.NE.1  .AND. MFR.NE.3  .AND. MFR.NE.5  .AND.
     &        MFR.NE.7  .AND. MFR.NE.9  .AND. MFR.NE.73 .AND.
     &        MFR.NE.27 .AND. MFR.NE.75 .AND. MFR.NE.79 ) THEN
            CALL ERREUR(16)
            GOTO 99
          ENDIF
        ENDIF
 
C Formulation 'ELECTROSTATIQUE' :
C Petite verification (a priori sans probleme)
        IF (LESFOR(1) .EQ. 'ELECTROSTATIQUE ') THEN
          IF (MFR .NE. 1) THEN
            CALL ERREUR(21)
            GOTO 99
          ENDIF
        ENDIF
 
C +--------------------------------------------------------------------+
C | Remplissage INFMOD et INFELE du IM-eme modele elementaire IMODEL   |
C +--------------------------------------------------------------------+
        IF (NGINT.NE.0.AND.MELE.NE.28) THEN
          CALL ERREUR(608)
          GOTO 990
        ENDIF
        INFMOD(1) = NGINT
 
C (fdp) Pour les elements JOI1 seulement, on stocke -1*ILIE dans INFMOD(9)
        IF (ILIE.NE.0) THEN
          IF (MELE.NE.265) THEN
            CALL ERREUR(19)
            GOTO 990
          ENDIF
          INFMOD(9) = -1 * ILIE
        ENDIF
*  AM  cas non-local
        IF (INLOC.NE.0) THEN
          INFMOD(13) = -1*INLOC
          INFMOD(14) = LULVIA
        ENDIF
 
C  Initialisation du infele et des segments d'integration
        infele(2) = NGINT
        infele(3) = NGMAS
        infele(4) = NGCON
        infele(6) = NGRIG
 
C  Cas particulier des relations de conformite pour les SURE
        IF (ITYP1.EQ.48) THEN
          imodel.infele( 1) = nefmod
          imodel.infele(14) = 48
          imodel.infele(13) = mfr2
C* Serait-il interessant de stocker mele=nepapa dans infele ?
c-dbg      write(ioimp,*) 'ITYP1 = SURE',ityp1,nefmod,mele,mfr,mfr2
        ENDIF
 
        CALL prquoi(imodel)
        if (ierr.ne.0) return
 
C +--------------------------------------------------------------------+
C | Initialisation des nomid (NOMS des composantes)                    |
C +--------------------------------------------------------------------+
 
C     cas particulier des relations de conformite pour les SURE
c     on recupere les noms de composantes 'DEPLACEM' et 'FORCES'
c     des elements parents (NEPAPA => QUA4 ou CUB8)
        if (ITYP1.EQ.48) then
          SEGINI,imode5=IMODEL
          imode5.NEFMOD=NEPAPA
          call inomid(imode5,lucvar,lucmat,lucmaf,luparx)
          if (ierr.ne.0) return
          imodel.LNOMID(1) = imode5.LNOMID(1)
          imodel.LNOMID(2) = imode5.LNOMID(2)
          SEGSUP,imode5
        else
          call inomid(imodel,lucvar,lucmat,lucmaf,luparx)
          if (ierr.ne.0) return
        endif
 
C       Test CLEMENT entre INFELE(16) et la dimension du NOMID des DEFORMATIONS
C       ATTENTION (celui des CONTRAINTES peut contenir une info en plus sur les MODES en fourier...)
        nomid = imodel.LNOMID(5)
        IF (nomid.GT.0) THEN
          imodel.INFELE(16) = nomid.LESOBL(/2) + nomid.LESFAC(/2)
        ELSE
          imodel.INFELE(16) = 0
        ENDIF
 
C +--------------------------------------------------------------------+
C | Quelques verifications supplementaires                             |
C +--------------------------------------------------------------------+
 
C=DEB==== FORMULATION HHO ==== Verification des noms primales/duales====
        IF (loHHO) THEN
          nomid1 = imodel.LNOMID(1)
          nomid2 = imodel.LNOMID(2)
c*          SEGACT,nomid1,nomid2
          n_z1 = nomid1.LESOBL(/2)
          n_z2 = nomid2.LESOBL(/2)
          IF (n_z1.EQ.0 .OR. n_z1.NE.n_z2) THEN
            write(ioimp,*) 'MODELI HHO: PRIMAL/DUAL number incorrect'
            CALL ERREUR(5)
            RETURN
          END IF
          DO i = 1, n_z1
            CALL VERMDI(nomid1.LESOBL(i),nomid2.LESOBL(i))
            IF (IERR.NE.0) RETURN
          END DO
          n_z1 = nomid1.LESFAC(/2)
          n_z2 = nomid2.LESFAC(/2)
          IF (n_z1.NE.0 .OR. n_z2.NE.0) THEN
            write(ioimp,*) 'MODELI HHO: LESFAC incorrect'
            CALL ERREUR(5)
            RETURN
          END IF
        END IF
C=FIN==== FORMULATION HHO ==============================================
 
C - Modele UTILISATEUR :
C   Verification que les composantes materiaux "obligatoires" sont declarees
        IF (LDIEXT) THEN
          NOMID  = IMODEL.LNOMID(6)
          NBROBL = NOMID.LESOBL(/2)
          MLMOT1 = lucmat
          SEGACT,MLMOT1
          NBCOMP = MLMOT1.MOTS(/2)
          ICOMP  = 0
          DO i = 1, NBROBL
            CALL PLACE(MLMOT1.MOTS,NBCOMP,IPLAC,NOMID.LESOBL(i))
            IF (IPLAC.EQ.0) THEN
              WRITE(IOIMP,80) MOTS(i)
 80           FORMAT('La composante obligatoire "',A8,'" est absente')
            ELSE
              ICOMP = ICOMP+1
            ENDIF
          ENDDO
          IF (ICOMP.NE.NBROBL) GOTO 99
        ENDIF
 
        mfr2 = INFELE(13)
        IF (FORMOD(1).EQ.'CONTRAINTE') mfr2 = 0
        ipmo = imodel
        CALL cotemo(ipmo,mfr2)
        IF (IERR.NE.0) RETURN
 
C +--------------------------------------------------------------------+
C | Point support pour les modes en defo. GENE (IFOUR=-3, 7 a 11, 14)
C | Ce point n'est pris en compte que si cela est necessaire
C +--------------------------------------------------------------------+
        mfr3 = mfr2
        CALL INFDPG(mfr3,IFOUR, LOGRE,ndpge)
        IF (LOGRE) THEN
C Erreur si ce point support n'est pas fourni avec le mot-cle GENE.
          IF (IPTGEN.EQ.0) THEN
            CALL ERREUR(925)
            RETURN
          ENDIF
          imodel.IPDPGE = IPTGEN
        ELSE
          IF (IPTGEN.NE.0) THEN
            write(ioimp,*) 'Mot-cle GENE + Point ignores...'
          ENDIF
          imodel.IPDPGE = 0
        ENDIF
 
 
        SEGACT,IMODEL*NOMOD
 
 10   CONTINUE
C ****************************************************************
C  Fin de la boucle (10) sur les maillages elementaires de IPGEOM
C ****************************************************************
 
* En cas de contact symetrique, on met tout dans le meme modele
      n1o = kmodel(/1)
      n1 = n1o
      do i = 1, n1o
        imode1 = mmode2.kmodel(i)
        if (imode1.ne.0) then
          n1 = n1+1
        endif
      enddo
* On a trouve du contact :
      if (n1.gt.n1o) then
        segadj mmodel
        nsou1 = n1
        do i = 1, n1o
          imode1 = mmode2.kmodel(i)
          if (imode1.ne.0) then
            kmodel(n1)=imode1
            n1=n1-1
            imodel=kmodel(i)
            imode1.nefmod=nefmod
            imode1.conmod=conmod
            do ip=1,infmod(/1)
              imode1.infmod(ip)=infmod(ip)
            enddo
            do ip=1,formod(/2)
              imode1.formod(ip)=formod(ip)
            enddo
            do ip=1,matmod(/2)
              imode1.matmod(ip)=matmod(ip)
            enddo
            imode1.ipdpge=ipdpge
            imode1.cmatee=cmatee
            imode1.imatee=imatee
            imode1.inatuu=inatuu
            imode1.ideriv=ideriv
            do ip=1,lnomid(/1)
              imode1.lnomid(ip)=lnomid(ip)
            enddo
            do ip=1,infele(/1)
              imode1.infele(ip)=infele(ip)
            enddo
            do ip=1,tymode(/2)
              imode1.tymode(ip)=tymode(ip)
            enddo
          endif
        enddo
        n1 = nsou1
      endif
      segsup mmode2
 
      IPMODE=MMODEL
 
C en cas de modele STAT ddddddddddddddddddddddddddddddddddddddddddd
  91  IF (IPTABS.GT.0.OR.IPMOD1.GT.0) THEN
c verification formulation
       IF (ipmod1.gt.0) THEN
         CALL ACTOBJ('MMODEL',IPMOD1,1)
         if (ierr.ne.0) return
         mmode1 = ipmod1
         imode1 = mmode1.kmodel(1)
         do jj=1,NFOR
         if (imode1.formod(jj).ne.LESFOR(JJ)) then
           call erreur(21)
           return
         endif   
         enddo
       ENDIF
c duplique le modele cree       
       if (ipmod1.le.0) ipmod1 = ipmode
C modele : pointer le modele elementaire approprie
       IF (iptabm.eq.0) THEN
         MMODE1 = ipmod1
         DO im = 1,kmodel(/1)
          imodel = kmodel(im)
          segact imodel*mod
          nobmod = ivamod(/1)
          nobmod = nobmod + 1
            nfor = formod(/2)
            nmat = matmod(/2)
            mn3 = infmod(/1)
          segadj imodel
          kbmod = 0
          do im1 = 1,MMODE1.KMODEL(/1)
            imode1 = mmode1.kmodel(im1)
            imomo = imode1
            lostat = .true.
 
C  criteres de verif assez sommaires ...
           if (imode1.nefmod.eq.nefmod.and.
     & imode1.imamod.ne.imamod.and.
     & (imode1.matmod(/2).eq.matmod(/2).or.
     &    imode1.matmod(/2).eq.(matmod(/2)-1)).and.
     & imode1.formod(/2).eq.formod(/2)) then
             do lmo = 1,formod(/2)
               if (formod(lmo).ne.imode1.formod(lmo)) lostat = .false.
             enddo
             do lmo = 1,imode1.matmod(/2)
              if (matmod(lmo).ne.imode1.matmod(lmo)) lostat = .false.
             enddo
           else
             lostat = .false.
           endif
           if (lostat.and.formod(1).eq.'MELANGE') then
C  verifs supplementaires : les modeles de ivamod sont ils bien construi
             lomela = .true.
             if ((nobmod - imode1.ivamod(/1)).gt.1) lomela = .false.
             if (imode1.ivamod(/1).gt.0) then
              do ivm3 = 1,imode1.ivamod(/1)
                  IF(imode1.tymode(ivm3).eq.'IMODEL') THEN
                    imode3 = imode1.ivamod(ivm3)
                    segact imode3
                  ENDIF
              enddo
             endif
             IF (nobmod.gt.1) THEN
             do ivm1 = 1,(nobmod-1)
               if (tymode(ivm1).eq.'IMODEL  ') then
              imode2 = ivamod(ivm1)
              segact imode2
cc
              if (imode2.ivamod(/1).ge.1) then
              do ivm2 = 1,imode2.ivamod(/1)
               if (imode2.tymode(ivm2).eq.'IMODEL') then
                 imode4 = imode2.ivamod(ivm2)
                 segact imode4
                 if (imode1.ivamod(/1).ge.1) then
                 do ivm3 = 1,imode1.ivamod(/1)
                  IF(imode1.tymode(ivm3).eq.'IMODEL') THEN
                    imode3 = imode1.ivamod(ivm3)
cc
             lostat = .true.
C  criteres de verif assez faibles ...
            if (imode3.nefmod.eq.imode4.nefmod.and.
     & imode3.imamod.eq.imode4.imamod.and.
     & imode3.matmod(/2).eq.imode4.matmod(/2).and.
     & imode3.conmod(17:24).eq.imode4.conmod(17:24).and.
     & imode3.formod(/2).eq.imode4.formod(/2)) then
             do lmo = 1,imode4.formod(/2)
       if (imode4.formod(lmo).ne.imode3.formod(lmo)) lostat = .false.
             enddo
             do lmo = 1,imode4.matmod(/2)
       if (imode4.matmod(lmo).ne.imode3.matmod(lmo)) lostat = .false.
             enddo
            else
             lostat = .false.
            endif
            if (lostat) then
              goto 75
            endif
cc
                  ENDIF
                 enddo
                 else
                  lostat = .false.
                 endif
               endif
              enddo
C
              else
               lomela = .false.
              endif
 
 75           lomela = lomela.and.lostat
               endif
 
             enddo
             ENDIF
             lostat = lomela
             do ivm3 = 1,imode1.ivamod(/1)
c               imode1 = imomo
                  IF(imode1.tymode(ivm3).eq.'IMODEL') THEN
                    imode3 = imode1.ivamod(ivm3)
                  ENDIF
             enddo
           endif
           if (lostat) then
             kbmod = kbmod + 1
             tymode(nobmod) = 'IMODEL'
             ivamod(nobmod) = imomo
             goto 79
           endif
          enddo
C *** ca se passe mal
          if (kbmod.ne.1) then
              write(ioimp,*) ' STATIO EN DEFAUT voir notice ',kbmod,im
              call erreur(251)
              goto 990
          endif
C ***
 79      CONTINUE
         ENDDO
       ENDIF
 
C : table : dupliquer modele elementaire et pointer
       if (iptabm.gt.0) then
        call modsta(ipmode,iptabm,ipmod1)
       endif
 
      ENDIF
C  fin du modele STAT ddddddddddddddddddddddddddddddddddddddddddddddd
 
      if (plicon.ne.0) segsup,plicon
 
C  Ecriture de l'objet MODELE cree
      CALL ACTOBJ('MMODEL  ',IPMODE,1)
      CALL ECROBJ('MMODEL  ',IPMODE)
      RETURN
 
C  Traitement des ERREURS
 99   CONTINUE
      CALL ERREUR(21)
 990  CONTINUE
      DO im = 1, kmodel(/1)
        imodel = kmodel(im)
        IF (imodel.NE.0) SEGSUP,imodel
      ENDDO
      SEGSUP,MMODEL
 
      if (plicon.ne.0) segsup,plicon
 
c      RETURN
      END