      SUBROUTINE ELIMIN
C
      INCLUDE 'com_coor.f'
      INCLUDE 'com_faces.f'
      INCLUDE 'com_options.f'
C
      INTEGER NODE(27),IPP(4,4,9),NUMF(6),IP8(4,6),NOD(8),IP827(8,8)
     &       ,IP9(4,4),IP4(2,4),IP4T(3,4),IP6(3,4)
     &       ,NFPR(6),IBID(28),NONO(8),IFAFA(6)
      REAL*4  PRES(6),XBID(28)
      EQUIVALENCE(IBID(2),NODE)
      EQUIVALENCE(IBID,XBID)
      EQUIVALENCE(NOD,NODE)
      LOGICAL*4 ZTRI,EGAL
      CHARACTER*128 CBID
      DATA IPP / 25,18, 6,13,  25,17, 1, 9,
     &           25, 9, 2,18,  25,13, 5,17,
     &           23,14, 7,15,  23,16, 5,13,
     &           23,13, 6,14,  23,15, 8,16,
     &           22,19, 7,14,  22,18, 2,10,
     &           22,10, 3,19,  22,14, 6,18,
     &           21,11, 3,10,  21, 9, 1,12,
     &           21,12, 4,11,  21,10, 2, 9,
     &           26,15, 7,19,  26,11, 4,20,
     &           26,20, 8,15,  26,19, 3,11,
     &           24,16, 8,20,  24,12, 1,17,
     &           24,17, 5,16,  24,20, 4,12,
     &           27,23,14,22,  27,24,16,23,
     &           27,21,12,24,  27,22,10,21,
     &           27,26,15,23,  27,23,13,25,
     &           27,25, 9,21,  27,21,11,26,
     &           27,22,19,26,  27,26,20,24,
     &           27,24,17,25,  27,25,18,22/
      DATA NUMF / 3,6,2,5,4,1 /
      DATA IP4 / 4,1 , 2,3 , 1,2 , 3,4 /
      DATA IP4T / 1,3,2
     &           ,2,3,4
     &           ,1,2,4
     &           ,1,4,3 /
      DATA IP6  / 6,1,4
     &           ,4,2,5
     &           ,5,3,6
     &           ,4,5,6 /
      DATA IP8  / 5,8,4,1
     &           ,2,3,7,6
     &           ,5,1,2,6
     &           ,4,8,7,3
     &           ,1,4,3,2
     &           ,8,5,6,7  /
      DATA IP9  / 9,5,2,6
     &           ,9,6,3,7
     &           ,9,7,4,8
     &           ,9,8,1,5  /
      DATA IP827 /  1, 9,21,12,17,25,27,24
     &            , 9, 2,10,21,25,18,22,27
     &            ,21,10, 3,11,27,22,19,26
     &            ,12,21,11, 4,24,27,26,20
     &            ,17,25,27,24, 5,13,23,16
     &            ,25,18,22,27,13, 6,14,23
     &            ,27,22,19,26,23,14, 7,15
     &            ,24,27,26,20,16,23,15, 8  /
      DATA US3 / .33333333333333333333 /
C
      IF (ISTDOUT.EQ.0) THEN
        IF (ILANG.EQ.0) THEN
          PRINT*,'Fin de lecture des noeuds. Lecture des éléments...'
        ELSE
          PRINT*,'End of nodes reading. Beginning elements reading..'
        ENDIF
      ENDIF
      DO I=1,NUMNP
        ISD(I) = 0
        INDEX(I) = 0
      ENDDO
      IF (NBCORN.NE.0) THEN
        DO I=1,NUMNP
          RAYON2(I) = 0.
        ENDDO
      ENDIF
C
      NDSMAX = 0
      IFACEXT = 0
      IVOL   = 0
      NEL    = 0
      ICPT   = 0
      NFC    = 0
      IGROUP = 0
      NUMSD  = 0
      NUMEL  = 0
      NUMNP2 = NUMNP
      MODTET = 0
      ICOIN  = 0
      IVIEUX = 0
 1000 CONTINUE
      IF (ICOURB.EQ.-5) THEN
        IGROUP = 1
        ICPT0 = 1
        NELEM = NELDEL
        NEL = NELEM
        IF (ICOURXYZ.EQ.0) THEN
          NDS = 3
          DO N=1,NEL
            NN = (N-1)*NDS+1
            CALL QUALITY(NODEL(NN),N,NDS)
          ENDDO
        ELSE
          NDS = 2
        ENDIF
        NDSEL = NDS
      ELSE
        IF (IARCH.EQ.0) THEN
          READ(IFAC,END=9999,ERR=9999) N,NTYP,NELEM,NDS
        ELSE
          CALL litrecbin(IFAC,NODE,LLL,0)
          IF (LLL.LE.0) GOTO 9999
          N     = NODE(1)
          NTYP  = NODE(2)
          NELEM = NODE(3)
          NDS   = NODE(4)
        ENDIF
        IF (NELEM.LE.0) GOTO 1000
        NDSEL = NDS
        IF (NDS*(NEL+NELEM).GT.NEMAX) CALL TROPDEPOINTS(NEL+NELEM,NDS,0)
CC maintenant on lit le no de sd directement
CCC        IGROUP = IGROUP+1
        IGROUP = N
        ICPT0  = ICPT+1
        NEL = NEL + NELEM
      ENDIF
      NMODP = MAX(1,NELEM/50)
      IF (NMODP.EQ.1) NMODP = NFMAX
C
C element a 27 noeuds
C
      IF (NDS.GE.27) THEN
        IF (ILANG.EQ.0) THEN
          ELEMENTS = 'Hexaèdres 27 noeuds'
        ELSE
          ELEMENTS = 'Hexaedrons 27 nodes'
        ENDIF
        MODTET = 40
        LELEM  = 19
        IOPREF = 1
        NBFAC  = 24
        NFEXT  = 6
        ICOIN  = 1
        IOPT   = 0
        DO N=1,NELEM
          IF (MOD(N,NMODP).EQ.0.AND.ISTDOUT.EQ.0)
     &         CALL MYFLUSH(100.*REAL(N)/REAL(NELEM))
          IF (IARCH.EQ.0) THEN
            READ(IFAC) NE,(NODE(J),J=1,NDS)
          ELSE
            CALL litrecbin(IFAC,IBID,LLL,0)
          ENDIF
          IMAT = IREF3(NODE(27))
          NUMEL = NUMEL+1
          CALL ECRNOD(NODE,NODEL,NUMEL,NDS)
          DO I=1,NDS
            IF (ISD(NODE(I)).EQ.0) THEN
              ISD(NODE(I)) = IGROUP
            ELSEIF(ISD(NODE(I)).NE.IGROUP) THEN
              ISD(NODE(I)) = -1
            ENDIF
          ENDDO
          DO I=9,20
            INDEX(NODE(I)) = INDEX(NODE(I))+1
          ENDDO
          DO I=21,26
            INDEX(NODE(I)) = -1
          ENDDO
C
C Decomposition de l'element a 27 noeuds en 24 faces exterieures
C
          DO IFCC = 1,NFEXT
            IF (IELIMI.EQ.0) IOPT = 1
            NUMFAC = NUMF(IFCC)
            DO K=1,4
              CALL TRIAGE(NODE(IPP(1,K,IFCC))
     &                   ,NODE(IPP(2,K,IFCC))
     &                   ,NODE(IPP(3,K,IFCC))
     &                   ,NODE(IPP(4,K,IFCC))
     &                   ,ICPT,ICPT0,NFC,IOPT,NUMEL,NUMFAC,NODE(27))
              ISD2(ICPT) = IGROUP
              ISD3(ICPT) = IMAT
            ENDDO
            IF (NFC.GT.NFCM) CALL COMPAC(ICPT,ICPT0,NFC)
          ENDDO
C
C Calcul des rayons de courbure et puissances dioptriques pour les yeux
C
          DO K=1,6
            IF (ICOR(NODE(20+K)).NE.0) CALL RAYFAC(NODE,K)
          ENDDO
C
C Decomposition en tetraedres pour isosurfaces
C
          DO K=1,8
            DO KK=1,8
              NONO(KK) = NODE(IP827(KK,K))
            ENDDO
            CALL TETRA(NONO)
          ENDDO 
        ENDDO
C
C 4 noeuds....
C
      ELSEIF(NDS.EQ.4) THEN
C
C ... 2d --> quadrangles
C
        IF (I2D.NE.0.OR.NTYP.EQ.3) THEN
          NBFAC = 1
          IPREFC = -1
          IOPREF = 0
          IF (I2D.EQ.0) THEN
            ELEMENTS = 'Quadrangles 3D Q1'
          ELSE
            ELEMENTS = 'Quadrangles 2D Q1'
          ENDIF
          LELEM = 17
          DO N=1,NELEM
          IF (MOD(N,NMODP).EQ.0.AND.ISTDOUT.EQ.0)
     &           CALL MYFLUSH(100.*REAL(N)/REAL(NELEM))
            ICPT = ICPT+1
            IF (IARCH.EQ.0) THEN
              IF (IVIEUX.EQ.0) THEN
                READ(IFAC,END=3131,ERR=3131) NE,(NOD(J),J=1,NDS)
     &                                    ,NPR,(NFPR(K),PRES(K),K=1,NPR)
                GOTO 3132
 3131           NPR = 0
                IVIEUX = 1
 3132           NUMEL = NUMEL+1
              ELSE
                READ(IFAC) NE,(NOD(J),J=1,NDS)
                NUMEL = NUMEL+1
              ENDIF
            ELSE
              CALL litrecbin(IFAC,IBID,LLL,0)
cc              NE = IBID(1)
              IF (LLL/4.GT.NDS+1) THEN
                NPR = IBID(NDS+2)
                IF (NPR.GT.0) THEN
                  DO J=1,NPR
                    NFPR(J) = IBID(NDS+1+2*J)
                    PRES(J) = XBID(NDS+2*(J+1))
                  ENDDO
                ENDIF
              ELSE
                NPR = 0
              ENDIF
              NUMEL = NUMEL+1
            ENDIF
            IF (NPR.GT.0) THEN
              DO J=1,NPR
                IF (PRES(J).EQ.9999.) THEN
                  IR = 5
                ELSEIF(PRES(J).EQ.-9999.) THEN
                  IR = 1
                ELSE
                  IR = MOD(NINT(PRES(J)),16)
                  IF (IR.LT.0) IR = IR+16
                ENDIF
                IREF(NOD(IP4(1,NFPR(J)))) = IR
                IREF(NOD(IP4(2,NFPR(J)))) = IR
              ENDDO
            ENDIF
            DO J=1,NDS
              ICLAS(J,ICPT) = NOD(J)
            ENDDO
            CALL ECRNOD(ICLAS(1,ICPT),NODEL,NUMEL,NDS)
            DO J=1,NDS
              NN = ICLAS(J,ICPT)
              INDEX(NN) = INDEX(NN)+1
              IF (INDEX(NN).LE.NVMX) NVOI(INDEX(NN),NN) = ICPT
              IF (ISD(NN).EQ.0) THEN
                ISD(NN) = IGROUP
              ELSEIF(ISD(NN).NE.IGROUP) THEN
                ISD(NN) = -1
              ENDIF
            ENDDO
            ICLAS(5,ICPT) = NUMEL
            ICLAS(6,ICPT) = 0
            NUMNP2 = NUMNP2+1
            X(NUMNP2) = .25*( X(ICLAS(1,ICPT))+X(ICLAS(2,ICPT))
     &                       +X(ICLAS(3,ICPT))+X(ICLAS(4,ICPT)) )
            Y(NUMNP2) = .25*( Y(ICLAS(1,ICPT))+Y(ICLAS(2,ICPT))
     &                       +Y(ICLAS(3,ICPT))+Y(ICLAS(4,ICPT)) )
            Z(NUMNP2) = .25*( Z(ICLAS(1,ICPT))+Z(ICLAS(2,ICPT))
     &                       +Z(ICLAS(3,ICPT))+Z(ICLAS(4,ICPT)) )
            DEPX(NUMNP2) = .25*(DEPX(ICLAS(1,ICPT))+DEPX(ICLAS(2,ICPT))
     &                         +DEPX(ICLAS(3,ICPT))+DEPX(ICLAS(4,ICPT)))
            DEPY(NUMNP2) = .25*(DEPY(ICLAS(1,ICPT))+DEPY(ICLAS(2,ICPT))
     &                         +DEPY(ICLAS(3,ICPT))+DEPY(ICLAS(4,ICPT)))
            DEPZ(NUMNP2) = .25*(DEPZ(ICLAS(1,ICPT))+DEPZ(ICLAS(2,ICPT))
     &                         +DEPZ(ICLAS(3,ICPT))+DEPZ(ICLAS(4,ICPT)))
            ICLAS(7,ICPT) = NUMNP2
            ISD2(ICPT)    = IGROUP
          ENDDO
C
C ... 3d --> tetraedres
C
        ELSE
          IOPREF = 1
          NBFAC = 4
          MODTET = 1
          IF (ILANG.EQ.0) THEN
            ELEMENTS = 'Tetraèdres'
            LELEM    = 10
          ELSE
            ELEMENTS = 'Tetraedrons'
            LELEM    = 11
          ENDIF
          IF (IELIMI.EQ.0) THEN
            IOPT = 4
          ELSE
            IOPT = 0
          ENDIF
          IF (NEL.EQ.NELEM) CALL CHECKFACEXT(CBID,LBID,IFACEXT)
          IBORLU = 0
 4436     CONTINUE
          DO N=1,NELEM
            IF (MOD(N,NMODP).EQ.0.AND.ISTDOUT.EQ.0)
     &           CALL MYFLUSH(100.*REAL(N)/REAL(NELEM))
            IF (IFACEXT.GT.0.AND.IBORLU.EQ.0) THEN
              READ(IFACEXT,*,END=4434,ERR=4434)
     &             NB,NFA,(IFAFA(I),I=1,NFA)
              IBORLU = 1
              GOTO 4435
 4434         IF (N.EQ.1) THEN
                IFACEXT = -1
                IF (ILANG.EQ.0) THEN
                  PRINT*,'*** Erreur à la lecture de '//CBID(1:LBID)//
     &                 ' - élément',N
                ELSE
                  PRINT*,'*** Error reading '//CBID(1:LBID)//
     &                 ' - element',N
                ENDIF
                NUMEL  = 0
                NUMNP2 = NUMNP
                GOTO 4436
              ENDIF
 4435         CONTINUE
            ENDIF
            IF (IARCH.EQ.0) THEN
              READ(IFAC) NE,N1,N2,N3,N4
            ELSE
              CALL litrecbin(IFAC,IBID,LLL,0)
cc              NE = IBID(1)
              N1 = IBID(2)
              N2 = IBID(3)
              N3 = IBID(4)
              N4 = IBID(5)
            ENDIF
            IVOL = IVOL+1
            NTET(1,IVOL) = N1
            NTET(2,IVOL) = N2
            NTET(3,IVOL) = N3
            NTET(4,IVOL) = N4
            NUMEL = NUMEL+1
            CALL ECRNOD(NTET(1,IVOL),NODEL,NUMEL,4)
            NUMNP2 = NUMNP2+1
            X(NUMNP2) = .25*( X(N1)+X(N2)+X(N3)+X(N4) )
            Y(NUMNP2) = .25*( Y(N1)+Y(N2)+Y(N3)+Y(N4) )
            Z(NUMNP2) = .25*( Z(N1)+Z(N2)+Z(N3)+Z(N4) )
            DEPX(NUMNP2) = .25*( DEPX(N1)+DEPX(N2)+DEPX(N3)+DEPX(N4) )
            DEPY(NUMNP2) = .25*( DEPY(N1)+DEPY(N2)+DEPY(N3)+DEPY(N4) )
            DEPZ(NUMNP2) = .25*( DEPZ(N1)+DEPZ(N2)+DEPZ(N3)+DEPZ(N4) )
            IF (IFACEXT.GT.0) THEN
              IF (NB.EQ.NUMEL) THEN
                IBORLU = 0
                DO K=1,4
                  IFLAG = 0
                  DO KK=1,NFA
                    IF (IFAFA(KK).EQ.K) IFLAG=KK
                  ENDDO
                  IF (IFLAG.NE.0) THEN
                    IOPT = 0
                    CALL ORDON2(NTET(IP4T(1,K),IVOL)
     &                   ,NTET(IP4T(2,K),IVOL)
     &                   ,NTET(IP4T(3,K),IVOL),M1,M2,M3)
                    CALL TRIAGE(M1,M2,M3,0,ICPT,ICPT0,NFC,IOPT,NUMEL,K
     &                   ,NUMNP2)
                    ISD2(ICPT) = IGROUP
                  ENDIF
                ENDDO
              ENDIF
            ELSE
              ICPT00 = ICPT
              CALL ORDON2(N1,N4,N3,M1,M2,M3)
              CALL TRIAGE(M1,M2,M3,0,ICPT,ICPT0,NFC,IOPT,NUMEL,4,NUMNP2)
              ISD2(ICPT) = IGROUP
              CALL ORDON2(N1,N2,N4,M1,M2,M3)
              CALL TRIAGE(M1,M2,M3,0,ICPT,ICPT0,NFC,IOPT,NUMEL,3,NUMNP2)
              ISD2(ICPT) = IGROUP
              CALL ORDON2(N2,N3,N4,M1,M2,M3)
              CALL TRIAGE(M1,M2,M3,0,ICPT,ICPT0,NFC,IOPT,NUMEL,2,NUMNP2)
              ISD2(ICPT) = IGROUP
              CALL ORDON2(N1,N3,N2,M1,M2,M3)
              CALL TRIAGE(M1,M2,M3,0,ICPT,ICPT0,NFC,IOPT,NUMEL,1,NUMNP2)
              ISD2(ICPT) = IGROUP
              DO K=ICPT00,ICPT
                DO J=1,NDS
                  NN = ICLAS(J,K)
                  INDEX(NN) = INDEX(NN)+1
                  IF (INDEX(NN).LE.NVMX) NVOI(INDEX(NN),NN) = K
                  JJ = ICLAS(J,K)
                  IF (ISD(JJ).EQ.0) THEN
                    ISD(JJ) = IGROUP
                  ELSEIF(ISD(JJ).NE.IGROUP) THEN
                    ISD(JJ) = -1
                  ENDIF
                ENDDO
              ENDDO
              IF (NFC.GT.NFCM) CALL COMPAC(ICPT,ICPT0,NFC)
            ENDIF
          ENDDO
        ENDIF
C
C Triangles P2
C
      ELSEIF(NDS.EQ.6) THEN
C
C c'est comme ca qu'on dit qu'on fait du Q2 (ou P2). idiot
C
        MODTET = 40
        IOPREF = 1
        NBFAC  = 4
        ICOIN  = 1
        IF (I2D.EQ.0) THEN
          ELEMENTS = 'Triangles 3D Q2'
        ELSE
          ELEMENTS = 'Triangles 2D Q2'
        ENDIF
        LELEM = 15
        DO N=1,NELEM
          IF (MOD(N,NMODP).EQ.0.AND.ISTDOUT.EQ.0)
     &         CALL MYFLUSH(100.*REAL(N)/REAL(NELEM))
          IF (IARCH.EQ.0) THEN
            READ(IFAC) NE,(NODE(J),J=1,NDS)
          ELSE
            CALL litrecbin(IFAC,IBID,LLL,0)
cc            NE = IBID(1)
          ENDIF
          NUMEL = NUMEL+1
          CALL ECRNOD(NODE,NODEL,NUMEL,NDS)
          DO I=1,4
            ICPT = ICPT+1
            DO J=1,3
              NN = NODE(IP6(J,I))
              INDEX(NN) = INDEX(NN)+1
              IF (INDEX(NN).LE.NVMX) NVOI(INDEX(NN),NN) = ICPT
              IF (ISD(NN).EQ.0) THEN
                ISD(NN) = IGROUP
              ELSEIF(ISD(NN).NE.IGROUP) THEN
                ISD(NN) = -1
              ENDIF
              ICLAS(J,ICPT) = NN
            ENDDO
            ICLAS(4,ICPT) = ICLAS(3,ICPT)
            ICLAS(5,ICPT) = NUMEL
            ICLAS(6,ICPT) = 0
            NUMNP2 = NUMNP2+1
            X(NUMNP2) = X(ICLAS(1,ICPT))
            Y(NUMNP2) = Y(ICLAS(1,ICPT))
            Z(NUMNP2) = Z(ICLAS(1,ICPT))
            DEPX(NUMNP2) = DEPX(ICLAS(1,ICPT))
            DEPY(NUMNP2) = DEPY(ICLAS(1,ICPT))
            DEPZ(NUMNP2) = DEPZ(ICLAS(1,ICPT))
            ICLAS(7,ICPT) = NUMNP2
            ISD2(ICPT)    = IGROUP
          ENDDO
        ENDDO
C
C Triangles P1 et segments
C
      ELSEIF(NDS.EQ.3.OR.NDS.EQ.2) THEN
ccc?        IOPREF = 0
        IOPREF = 1
        NBFAC = 1
        IF (NDS.EQ.3) THEN
          IF (I2D.EQ.0) THEN
            ELEMENTS = 'Triangles 3D P1'
          ELSE
            ELEMENTS = 'Triangles 2D P1'
          ENDIF
          LELEM = 15
        ELSE
          IF (I2D.EQ.0) THEN
            ELEMENTS = 'Segments 3D'
          ELSE
            ELEMENTS = 'Segments 2D'
          ENDIF
          LELEM = 11
        ENDIF
        DO N=1,NELEM
          IF (MOD(N,NMODP).EQ.0.AND.ISTDOUT.EQ.0)
     &         CALL MYFLUSH(100.*REAL(N)/REAL(NELEM))
          ICPT = ICPT+1
          IF (ICOURB.NE.-5) THEN
            IF (IARCH.EQ.0) THEN
              READ(IFAC) NE,(ICLAS(J,ICPT),J=1,NDS)
            ELSE
              CALL litrecbin(IFAC,IBID,LLL,0)
cc              NE = IBID(1)
              DO J=1,NDS
                ICLAS(J,ICPT) = IBID(J+1)
              ENDDO
            ENDIF
          ENDIF
          NUMEL = NUMEL+1
          CALL ECRNOD(ICLAS(1,ICPT),NODEL,NUMEL,NDS)
          ICLAS(4,ICPT) = ICLAS(3,ICPT)
          ICLAS(5,ICPT) = NUMEL
          ICLAS(6,ICPT) = 0
          NUMNP2 = NUMNP2+1
          IF (NDS.EQ.3) THEN
            X(NUMNP2) = ( X(ICLAS(1,ICPT))+X(ICLAS(2,ICPT))
     &                  + X(ICLAS(3,ICPT)) )*US3
            Y(NUMNP2) = ( Y(ICLAS(1,ICPT))+Y(ICLAS(2,ICPT))
     &                  + Y(ICLAS(3,ICPT)) )*US3
            Z(NUMNP2) = ( Z(ICLAS(1,ICPT))+Z(ICLAS(2,ICPT))
     &                  + Z(ICLAS(3,ICPT)) )*US3
            DEPX(NUMNP2) = ( DEPX(ICLAS(1,ICPT))+DEPX(ICLAS(2,ICPT))
     &                     + DEPX(ICLAS(3,ICPT)) )*US3
            DEPY(NUMNP2) = ( DEPY(ICLAS(1,ICPT))+DEPY(ICLAS(2,ICPT))
     &                     + DEPY(ICLAS(3,ICPT)) )*US3
            DEPZ(NUMNP2) = ( DEPZ(ICLAS(1,ICPT))+DEPZ(ICLAS(2,ICPT))
     &                     + DEPZ(ICLAS(3,ICPT)) )*US3
          ELSE
            X(NUMNP2) = ( X(ICLAS(1,ICPT))+X(ICLAS(2,ICPT)) )*.5
            Y(NUMNP2) = ( Y(ICLAS(1,ICPT))+Y(ICLAS(2,ICPT)) )*.5
            Z(NUMNP2) = ( Z(ICLAS(1,ICPT))+Z(ICLAS(2,ICPT)) )*.5
            DEPX(NUMNP2) = (DEPX(ICLAS(1,ICPT))+DEPX(ICLAS(2,ICPT)))*.5
            DEPY(NUMNP2) = (DEPY(ICLAS(1,ICPT))+DEPY(ICLAS(2,ICPT)))*.5
            DEPZ(NUMNP2) = (DEPZ(ICLAS(1,ICPT))+DEPZ(ICLAS(2,ICPT)))*.5
            ICLAS(3,ICPT) = ICLAS(2,ICPT)
          ENDIF
          ICLAS(7,ICPT) = NUMNP2
          ISD2(ICPT) = IGROUP
          DO J=1,NDS
            NN = ICLAS(J,ICPT)
            INDEX(NN) = INDEX(NN)+1
            IF (INDEX(NN).LE.NVMX) NVOI(INDEX(NN),NN) = ICPT
            JJ = ICLAS(J,ICPT)
            IF (ISD(JJ).EQ.0) THEN
              ISD(JJ) = IGROUP
            ELSEIF(ISD(JJ).NE.IGROUP) THEN
              ISD(JJ) = -1
            ENDIF
          ENDDO
        ENDDO
C
C Briques a 8 noeuds
C
      ELSEIF(NDS.EQ.8) THEN
        IOPREF = 0
        NBFAC  = 6
        MODTET = 5
        IF (ILANG.EQ.0) THEN
          ELEMENTS = 'Hexaèdres 8 noeuds'
        ELSE
          ELEMENTS = 'Hexaedrons 8 nodes'
        ENDIF
        LELEM    = 18
        IF (IELIMI.EQ.0) THEN
          IOPT = 8
        ELSE
          IOPT = 0
        ENDIF
        IF (NEL.EQ.NELEM) CALL CHECKFACEXT(CBID,LBID,IFACEXT)
        IBORLU = 0
 3436   CONTINUE
        DO N=1,NELEM
          IF (MOD(N,NMODP).EQ.0.AND.ISTDOUT.EQ.0)
     &         CALL MYFLUSH(100.*REAL(N)/REAL(NELEM))
          IF (IFACEXT.GT.0.AND.IBORLU.EQ.0) THEN
            READ(IFACEXT,*,END=3434,ERR=3434) NB,NFA,(IFAFA(I),I=1,NFA)
            IBORLU = 1
            GOTO 3435
 3434       IF (N.EQ.1) THEN
              IFACEXT = -1
              IF (ILANG.EQ.0) THEN
                PRINT*,'*** Erreur à la lecture de '//CBID(1:LBID)//
     &               ' - élément',N
              ELSE
                PRINT*,'*** Error reading '//CBID(1:LBID)//
     &               ' - element',N
              ENDIF
              NUMEL  = 0
              NUMNP2 = NUMNP
              GOTO 3436
            ENDIF
 3435       CONTINUE
          ENDIF
          IF (IARCH.EQ.0) THEN
            IF (IVIEUX.EQ.0) THEN
              READ(IFAC,END=3232,ERR=3232) NE,(NODE(J),J=1,NDS)
     &                                   ,NPR,(NFPR(J),PRES(J),J=1,NPR)
              GOTO 3233
 3232         NPR = 0
              IVIEUX = 1
 3233         NUMEL = NUMEL+1
            ELSE
              READ(IFAC) NE,(NODE(J),J=1,NDS)
              NUMEL = NUMEL+1
            ENDIF
          ELSE
            CALL litrecbin(IFAC,IBID,LLL,0)
cc            NE = IBID(1)
            IF (LLL/4.GT.NDS+1) THEN
              NPR = IBID(NDS+2)
              IF (NPR.GT.0) THEN
                DO J=1,NPR
                  NFPR(J) = IBID(NDS+1+2*J)
                  PRES(J) = XBID(NDS+2*(J+1))
                ENDDO
              ENDIF
            ELSE
              NPR = 0
            ENDIF
            NUMEL = NUMEL+1
          ENDIF
          IF (NPR.GT.0) THEN
            DO J=1,NPR
              IF (PRES(J).EQ.9999.) THEN
                IR = 5
              ELSEIF(PRES(J).EQ.-9999.) THEN
                IR = 1
              ELSE
                IR = MOD(NINT(PRES(J)),16)
                IF (IR.LT.0) IR = IR+16
              ENDIF
              IREF(NODE(IP8(1,NFPR(J)))) = IR
              IREF(NODE(IP8(2,NFPR(J)))) = IR
            ENDDO
          ENDIF
          CALL ECRNOD(NODE,NODEL,NUMEL,NDS)
          DO I=1,NDS
            IF (ISD(NODE(I)).EQ.0) THEN
              ISD(NODE(I)) = IGROUP
            ELSEIF(ISD(NODE(I)).NE.IGROUP) THEN
              ISD(NODE(I)) = -1
            ENDIF
          ENDDO
          NUMNP2 = NUMNP2+1
          X(NUMNP2) = 
     &      .125*( X(NODE(1))+X(NODE(2))+X(NODE(3))+X(NODE(4))
     &           + X(NODE(5))+X(NODE(6))+X(NODE(7))+X(NODE(8)) )
          Y(NUMNP2) = 
     &      .125*( Y(NODE(1))+Y(NODE(2))+Y(NODE(3))+Y(NODE(4))
     &           + Y(NODE(5))+Y(NODE(6))+Y(NODE(7))+Y(NODE(8)) )
          Z(NUMNP2) = 
     &      .125*( Z(NODE(1))+Z(NODE(2))+Z(NODE(3))+Z(NODE(4))
     &           + Z(NODE(5))+Z(NODE(6))+Z(NODE(7))+Z(NODE(8)) )
          DEPX(NUMNP2) = .125*( DEPX(NODE(1))+DEPX(NODE(2))
     &                        + DEPX(NODE(3))+DEPX(NODE(4))
     &                        + DEPX(NODE(5))+DEPX(NODE(6))
     &                        + DEPX(NODE(7))+DEPX(NODE(8)) )
          DEPY(NUMNP2) = .125*( DEPY(NODE(1))+DEPY(NODE(2))
     &                        + DEPY(NODE(3))+DEPY(NODE(4))
     &                        + DEPY(NODE(5))+DEPY(NODE(6))
     &                        + DEPY(NODE(7))+DEPY(NODE(8)) )
          DEPZ(NUMNP2) = .125*( DEPZ(NODE(1))+DEPZ(NODE(2))
     &                        + DEPZ(NODE(3))+DEPZ(NODE(4))
     &                        + DEPZ(NODE(5))+DEPZ(NODE(6))
     &                        + DEPZ(NODE(7))+DEPZ(NODE(8)) )
          IF (IFACEXT.GT.0) THEN
            IF (NB.EQ.NUMEL) THEN
              IBORLU = 0
              DO K=1,6
                IFLAG = 0
                DO KK=1,NFA
                  IF (IFAFA(KK).EQ.K) IFLAG=KK
                ENDDO
                IF (IFLAG.NE.0) THEN
                  IOPT = 0
                  CALL ORDONN(NODE(IP8(1,K)),NODE(IP8(2,K))
     &                       ,NODE(IP8(3,K)),NODE(IP8(4,K)),M1,M2,M3,M4)
                  CALL TRIAGE(M1,M2,M3,M4,ICPT,ICPT0,NFC,IOPT,NUMEL,K
     &                       ,NUMNP2)
                  ISD2(ICPT) = IGROUP
                ENDIF
              ENDDO
            ENDIF
          ELSE
            DO K=1,6
              CALL ORDONN(NODE(IP8(1,K)),NODE(IP8(2,K))
     &                   ,NODE(IP8(3,K)),NODE(IP8(4,K)),M1,M2,M3,M4)
              CALL TRIAGE(M1,M2,M3,M4,ICPT,ICPT0,NFC,IOPT,NUMEL,K
     &                   ,NUMNP2)
              ISD2(ICPT) = IGROUP
            ENDDO
            IF (NFC.GT.NFCM) CALL COMPAC(ICPT,ICPT0,NFC)
          ENDIF
          DO K=1,8
            INDEX(NODE(K)) = INDEX(NODE(K))+1
          ENDDO
          CALL TETRA(NODE)
        ENDDO
C
C 9 noeuds....
C
      ELSEIF(NDS.EQ.9) THEN
C
C c'est comme ca qu'on dit qu'on fait du Q2 (ou P2). idiot
C
        MODTET = 40
        IOPREF = 1
        NBFAC  = 4
        ICOIN  = 1
        IF (I2D.EQ.0) THEN
          ELEMENTS = 'Quadrangles 3D Q2'
        ELSE
          ELEMENTS = 'Quadrangles 2D Q2'
        ENDIF
        LELEM = 17
        DO N=1,NELEM
          IF (MOD(N,NMODP).EQ.0.AND.ISTDOUT.EQ.0)
     &         CALL MYFLUSH(100.*REAL(N)/REAL(NELEM))
          IF (IARCH.EQ.0) THEN
            READ(IFAC) NE,(NODE(J),J=1,NDS)
          ELSE
            CALL litrecbin(IFAC,IBID,LLL,0)
cc            NE = IBID(1)
          ENDIF
          IMAT = IREF3(NODE(9))
          NUMEL = NUMEL+1
          CALL ECRNOD(NODE,NODEL,NUMEL,NDS)
          DO I=1,4
            ICPT = ICPT+1
            DO J=1,4
              NN = NODE(IP9(J,I))
              INDEX(NN) = INDEX(NN)+1
              IF (INDEX(NN).LE.NVMX) NVOI(INDEX(NN),NN) = ICPT
              IF (ISD(NN).EQ.0) THEN
                ISD(NN) = IGROUP
              ELSEIF(ISD(NN).NE.IGROUP) THEN
                ISD(NN) = -1
              ENDIF
              ICLAS(J,ICPT) = NN
            ENDDO
            ICLAS(5,ICPT) = NUMEL
            ICLAS(6,ICPT) = 0
            NUMNP2 = NUMNP2+1
            X(NUMNP2) = X(ICLAS(1,ICPT))
            Y(NUMNP2) = Y(ICLAS(1,ICPT))
            Z(NUMNP2) = Z(ICLAS(1,ICPT))
            DEPX(NUMNP2) = DEPX(ICLAS(1,ICPT))
            DEPY(NUMNP2) = DEPY(ICLAS(1,ICPT))
            DEPZ(NUMNP2) = DEPZ(ICLAS(1,ICPT))
            ICLAS(7,ICPT) = NUMNP2
            ISD2(ICPT)    = IGROUP
            ISD3(ICPT)    = IMAT
          ENDDO
        ENDDO
      ELSE
        IF (ILANG.EQ.0) THEN
          PRINT*,'***** Eléments à',NDS,' noeuds de type inconnu *****'
        ELSE
          PRINT*,'*****',NDS,'-nodes elements are unsupported *****'
        ENDIF
        STOP
      ENDIF
      IF (ISTDOUT.EQ.0) CALL ENDFLUSH
      NDSMAX = MAX(NDSMAX,NDS)
      IF (ICOURB.NE.-5) GOTO 1000
C
C Recopie les coordonnees des points des faces selectionnees
C
 9999 IF (ICOURB.NE.-5) CLOSE(IFAC)
      IF (IFACEXT.GT.0) CLOSE(IFACEXT)
      ZTRI = NDSMAX.EQ.2.OR.NDSMAX.EQ.3.OR.NDSMAX.EQ.6
     &  .OR.(NDSMAX.EQ.4.AND.I2D.EQ.0.AND.NTYP.NE.3)
      IF (ZTRI) THEN
        NDSB = 3
      ELSE
        NDSB = 4
      ENDIF
      IF (INOPO.NE.0) CALL EXEC('/bin/rm -f '//NOM_FICH(1:LONG)//'*')
      NFACE = 0
      DO 30 N=1,ICPT
        IF (ICLAS(1,N).NE.0) THEN
          NFACE = NFACE + 1
          ISD2(NFACE) = ISD2(N)
          ISD3(NFACE) = ISD3(N)
          NUMSD = MAX(NUMSD,ISD2(NFACE))
          NNUMFA(NFACE) = ICLAS(5,N)
          NRFAC(NFACE)  = ICLAS(6,N) 
          DO I=1,NDSB
            NFAC(I,NFACE) = ICLAS(I,N)
            XF(I,NFACE) = X(ICLAS(I,N))
            YF(I,NFACE) = Y(ICLAS(I,N))
            ZF(I,NFACE) = Z(ICLAS(I,N))
          ENDDO
          NFAC(NDSB+1,NFACE) = ICLAS(7,N)
          XF(NDSB+1,NFACE) = X(ICLAS(7,N))
          YF(NDSB+1,NFACE) = Y(ICLAS(7,N))
          ZF(NDSB+1,NFACE) = Z(ICLAS(7,N))
          IF (ZTRI) THEN
            XMA = MAX(X(ICLAS(1,N)),X(ICLAS(2,N))
     &               ,X(ICLAS(3,N)))
            XMI = MIN(X(ICLAS(1,N)),X(ICLAS(2,N))
     &               ,X(ICLAS(3,N)))
            ISMIN = MIN( ISD(ICLAS(1,N)),ISD(ICLAS(2,N))
     &                  ,ISD(ICLAS(3,N)) )
            ISMAX = MAX( ISD(ICLAS(1,N)),ISD(ICLAS(2,N))
     &                  ,ISD(ICLAS(3,N)) )
          ELSE
            XMA = MAX(X(ICLAS(1,N)),X(ICLAS(2,N))
     &               ,X(ICLAS(3,N)),X(ICLAS(4,N)))
            XMI = MIN(X(ICLAS(1,N)),X(ICLAS(2,N))
     &               ,X(ICLAS(3,N)),X(ICLAS(4,N)))
            ISMIN = MIN( ISD(ICLAS(1,N)),ISD(ICLAS(2,N))
     &                  ,ISD(ICLAS(3,N)),ISD(ICLAS(4,N)) )
            ISMAX = MAX( ISD(ICLAS(1,N)),ISD(ICLAS(2,N))
     &                  ,ISD(ICLAS(3,N)),ISD(ICLAS(4,N)) )
          ENDIF
          IF (EGAL(XMI,XMA)) THEN
            IF (EGAL(XMI,-XMED).AND.ISYM.EQ.4) THEN
              IFPLAN(NFACE) = 3
            ELSE
              IFPLAN(NFACE) = 0
            ENDIF
          ELSE
            IF (ZTRI) THEN
              YMA = MAX(Y(ICLAS(1,N)),Y(ICLAS(2,N))
     &                 ,Y(ICLAS(3,N)))
              YMI = MIN(Y(ICLAS(1,N)),Y(ICLAS(2,N))
     &                 ,Y(ICLAS(3,N)))
            ELSE
              YMA = MAX(Y(ICLAS(1,N)),Y(ICLAS(2,N))
     &                 ,Y(ICLAS(3,N)),Y(ICLAS(4,N)))
              YMI = MIN(Y(ICLAS(1,N)),Y(ICLAS(2,N))
     &                 ,Y(ICLAS(3,N)),Y(ICLAS(4,N)))
            ENDIF
            IF (EGAL(YMI,YMA)) THEN
              IF (EGAL(YMI,-YMED)) THEN
                IFPLAN(NFACE) = 2
              ELSE
                IFPLAN(NFACE) = 0
              ENDIF
            ENDIF
          ENDIF
          IF (ISMIN.NE.ISMAX) THEN
            ITOTO = ISD2(NFACE)
            IP = 1000
            IF (ISD(ICLAS(1,N)).LT.0.AND.ISD(ICLAS(2,N)).LT.0) THEN
              ITOTO = ITOTO + IP
              IP = IP*10
            ENDIF
            IF (ISD(ICLAS(2,N)).LT.0.AND.ISD(ICLAS(3,N)).LT.0) THEN
              ITOTO = ITOTO + IP*2
              IP = IP*10
            ENDIF
            IF (ZTRI) THEN
              IF (ISD(ICLAS(3,N)).LT.0.AND.ISD(ICLAS(1,N)).LT.0) THEN
                ITOTO = ITOTO + IP*3
              ENDIF
            ELSE
              IF (ISD(ICLAS(3,N)).LT.0.AND.ISD(ICLAS(4,N)).LT.0) THEN
                ITOTO = ITOTO + IP*3
                IP = IP*10
              ENDIF
              IF (ISD(ICLAS(4,N)).LT.0.AND.ISD(ICLAS(1,N)).LT.0) THEN
                ITOTO = ITOTO + IP*4
              ENDIF
            ENDIF
            ISD2(NFACE) = ITOTO
          ENDIF
        ENDIF
 30   CONTINUE
      IGROUP = NUMSD
C
C Pourquoi c'etait commente ? 6/98
C
      IF (I2D.NE.0) NDS = NDSB
C
      IF (NDS.EQ.3.OR.NDS.EQ.2) IPREFC=-1
      IF (NDS.EQ.2) IFC = 1
      IF (ISTDOUT.EQ.0) THEN
        IF (ILANG.EQ.0) THEN
          IF (IVIEUX.NE.0)
     &       PRINT*,'*** Attention !! Vieux format de maillage.........'
          PRINT*,'Fin de lecture. Recherche de la frontière.........'
        ELSE
          IF (IVIEUX.NE.0)
     &       PRINT*,'*** Warning !! Old-style mesh file................'
          PRINT*,'End of reading. Computing the boundary............'
        ENDIF
      ENDIF
      CALL MINDIAM(NFACE)
      CALL FROETC(NTYP)
      NUMBIS = NUMNP2
      IF (ZTRI) NDS=3
C
C Determination des elements du bord
C
      DO I=1,NEL
        IELBOR(I) = 0
      ENDDO
      DO I=1,NFACE
        IF (IELBOR(NNUMFA(I)).EQ.0) THEN
          IELBOR(NNUMFA(I)) = I
        ELSE
          IELBOR(NNUMFA(I)) = -I
        ENDIF        
      ENDDO
      NNN = 0
      DO I=1,NEL
        IF (IELBOR(I).NE.0) NNN = NNN+1
      ENDDO
      IF (I2D.EQ.0.AND.IPS2D.EQ.0) THEN
        IF (ISTDOUT.EQ.0) THEN
          IF (ILANG.EQ.0) THEN
            PRINT*,NNN,' éléments de bord sur'
     &       ,NEL,' (',REAL(NINT(1000.*REAL(NNN)/REAL(NEL)))*.1,' %)'
          ELSE
            PRINT*,NNN,' boundary elements among'
     &       ,NEL,' (',REAL(NINT(1000.*REAL(NNN)/REAL(NEL)))*.1,' %)'
          ENDIF
        ENDIF
      ENDIF
C
C Ecriture du fichier .facext 
C (8 noeuds et 4 noeuds tetra seulement pour l'instant)
C
      IF (IFACEXT.EQ.-1.AND.IELIMI.EQ.0) THEN
        CALL PREMIER_LIBRE(IFACEXT)
        OPEN(IFACEXT,FILE=CBID(1:LBID),IOSTAT=IERR)
        IF (IERR.EQ.0) THEN
          DO I=1,NEL
            ITAB(I) = 0
          ENDDO
          DO I=1,NFACE
            ITAB(NNUMFA(I)) = ITAB(NNUMFA(I))+1
          ENDDO
          NN = 0
          DO I=1,NEL
            IF (ITAB(I).GT.0) THEN
              DO J=1,ITAB(I)
                NN = NN+1
                IFAFA(J) = NRFAC(NN)
              ENDDO
              WRITE(IFACEXT,*) I,ITAB(I),(IFAFA(K),K=1,ITAB(I))
            ENDIF
          ENDDO
          CLOSE(IFACEXT)
        ENDIF
      ENDIF
C
      END
C=======================================================================
      SUBROUTINE ECRNOD(NODE,NODEL,NUMEL,NDS)
C
      INTEGER NODE(*),NODEL(*)
C
      II = (NUMEL-1)*NDS
      DO I=1,NDS
        NODEL(II+I) = NODE(I)
      ENDDO
C
C calcul des criteres de qualite des EF
C
      CALL QUALITY(NODE,NUMEL,NDS)
      END
C=======================================================================
      SUBROUTINE QUALITY(NODE,N,ND)
      INCLUDE 'com_coor.f'
      INCLUDE 'com_faces.f'
C
      DIMENSION ANG(4)
      INTEGER   NODE(*),ISUC3(3),ISUC4(4),IPRE4(4),IP8(4,6)
      DATA ISUC3 / 2,3,1 /
      DATA ISUC4 / 2,3,4,1 /
      DATA IPRE4 / 4,1,2,3 /
      DATA IP8   / 5,8,4,1,
     &             2,3,7,6,
     &             5,1,2,6,
     &             4,8,7,3,
     &             1,4,3,2,
     &             8,5,6,7  /
C
      RMAX = 0.
      IF (ND.EQ.3.OR.ND.EQ.6) THEN
C
C Triangles
C
        XC = (X(NODE(1))+X(NODE(2))+X(NODE(3)))/3.
        YC = (Y(NODE(1))+Y(NODE(2))+Y(NODE(3)))/3.
        ZC = (Z(NODE(1))+Z(NODE(2))+Z(NODE(3)))/3.
        DO J=1,3
          K = ISUC3(J)
          L = ISUC3(K)
          ANG(J) = CALCANG(X,Y,Z,NODE,J,K,L)
          RR = SQRT((X(NODE(J))-XC)**2
     &             +(Y(NODE(J))-YC)**2
     &             +(Z(NODE(J))-ZC)**2)
          RMAX = AMAX1(RMAX,RR)
        ENDDO
        QUALITE(N,1) = AMIN1(ANG(1),ANG(2),ANG(3))
     &                /AMAX1(ANG(1),ANG(2),ANG(3))
      ELSEIF(ND.EQ.4.OR.ND.EQ.9) THEN
C
C quadrangles (4,9 noeuds) (a revoir pour les tetra)
C
        XC = (X(NODE(1))+X(NODE(2))+X(NODE(3))+X(NODE(4)))*0.25
        YC = (Y(NODE(1))+Y(NODE(2))+Y(NODE(3))+Y(NODE(4)))*0.25
        ZC = (Z(NODE(1))+Z(NODE(2))+Z(NODE(3))+Z(NODE(4)))*0.25
        DO J=1,4
          K = ISUC4(J)
          L = IPRE4(J)
          ANG(J) = CALCANG(X,Y,Z,NODE,J,K,L)
          RR = SQRT((X(NODE(J))-XC)**2
     &             +(Y(NODE(J))-YC)**2
     &             +(Z(NODE(J))-ZC)**2)
          RMAX = AMAX1(RMAX,RR)
        ENDDO
        QUALITE(N,1) = AMIN1(ANG(1),ANG(2),ANG(3),ANG(4))
     &                /AMAX1(ANG(1),ANG(2),ANG(3),ANG(4))
      ELSE
C
C Hexa 8, 27 noeuds
C
        ANGMIN = 5.
        ANGMAX = 0.
        XC = (X(NODE(1))+X(NODE(2))+X(NODE(3))+X(NODE(4))
     &       +X(NODE(5))+X(NODE(6))+X(NODE(7))+X(NODE(8)))*0.125
        YC = (Y(NODE(1))+Y(NODE(2))+Y(NODE(3))+Y(NODE(4))
     &       +Y(NODE(5))+Y(NODE(6))+Y(NODE(7))+Y(NODE(8)))*0.125
        ZC = (Z(NODE(1))+Z(NODE(2))+Z(NODE(3))+Z(NODE(4))
     &       +Z(NODE(5))+Z(NODE(6))+Z(NODE(7))+Z(NODE(8)))*0.125
        DO I=1,6
          DO II=1,4
            J = IP8(II,I)
            K = IP8(ISUC4(II),I)
            L = IP8(IPRE4(II),I)
            AA = CALCANG(X,Y,Z,NODE,J,K,L)
            ANGMAX = AMAX1(ANGMAX,AA)
            ANGMIN = AMIN1(ANGMIN,AA)
          ENDDO
        ENDDO
        QUALITE(N,1) = ANGMIN/ANGMAX
        DO J=1,8
          RR = SQRT((X(NODE(J))-XC)**2
     &             +(Y(NODE(J))-YC)**2
     &             +(Z(NODE(J))-ZC)**2)
          RMAX = AMAX1(RMAX,RR)
        ENDDO
      ENDIF
C
C Pas fini...
C
      QUALITE(N,2) = RMAX
C
      END
C=======================================================================
      REAL*4 FUNCTION CALCANG(X,Y,Z,NODE,J,K,L)
      INTEGER   NODE(*)
      DIMENSION X(*),Y(*),Z(*)
CC      DATA CQVSPI / 57.29577951 /
      DATA PI / 3.14159265358979  /
C
      JJ = NODE(J)
      KK = NODE(K)
      LL = NODE(L)
      U1 = X(KK)-X(JJ)
      V1 = Y(KK)-Y(JJ)
      W1 = Z(KK)-Z(JJ)
      U2 = X(LL)-X(JJ)
      V2 = Y(LL)-Y(JJ)
      W2 = Z(LL)-Z(JJ)
      XN1 = SQRT(U1**2+V1**2+W1**2)
      XN2 = SQRT(U2**2+V2**2+W2**2)
      IF (XN1*XN2.NE.0.) THEN
        ANG = ASIN(MIN(1.,SQRT((U1*V2-U2*V1)**2
     &                 +(V1*W2-V2*W1)**2
     &                 +(W1*U2-W2*U1)**2)/(XN1*XN2)))
        SCAL = U1*U2 + V1*V2 + W1*W2
        IF (SCAL.GE.0.) THEN
          CALCANG = ANG
        ELSE
          CALCANG = PI-ANG
        ENDIF
      ELSE
        CALCANG = -1.
      ENDIF
      END
