      SUBROUTINE CALSUR(IMP)
C
      INCLUDE 'com_coor.f'
      INCLUDE 'com_faces.f'
      INCLUDE 'com_options.f'
C
      REAL*8 S0,S1
      PARAMETER (NUMEPS=9)
      INTEGER IAR(2,6),IFA(3,4),ISUC(4),KCOMPL(4),ITMP(2),IG5(5)
      DIMENSION XTMP(2),YTMP(2),ZTMP(2),VTMP(2),EPSBES(NUMEPS)
      DATA IAR  / 1,2 , 2,3 , 3,1 , 1,4 , 2,4 , 3,4 /
      DATA IFA  / 1,3,2 , 1,2,4 , 2,3,4 , 1,4,3 /
      DATA KCOMPL / 4,3,1,2 /
      DATA ISUC / 2,3,4,1 /
      DATA EPSBES / 1.E-5 , 3.E-5 , 6.E-5 , 1.E-4 , 1.5E-4
     &             ,3.E-4 , 6.E-4 , 1.E-3 , 3.E-3 /
C
      ITOUS = 0
      DO I=1,NUMSD
        IF (ISDVU(I).EQ.0) ITOUS = ITOUS+1
      ENDDO
      EPST = 1.E-6
c      EPST = 1.E-5
c      EPST = 1.E-4
      UMEPST = 1.-EPST 
      IF (ISOCOUP.EQ.1) THEN
        EPSISO = EPST*(VMAXXYZ(IPCOUP)-VMINXYZ(IPCOUP))
      ELSE
        EPSISO = EPST*(VMAXISO-VMINISO)
      ENDIF
      visomi = VISO-EPSISO
      visoma = VISO+EPSISO
      CALL INV3X3(ROTA,ROTLOC,IERR)
      DO I=1,3
        DO J=1,3
          ROTA(J,I) = ROTLOC(J,I)
        ENDDO
      ENDDO
      CALL ROTATE(0)
      IPASS = 0
Cfj      NBIZB = IVOL
Cfj      IBEST = 0
Cfj      IDERN = 0
      IF (IMP.GT.0) CALL TEMPS(S0,IS)
C
    1 IPASS = IPASS+1
cc      EPS = EPSBES(IPASS)
      EPS = EPSISO
      EPS2 = EPS**2
C
    2 NBIZ0 = 0
      NBIZ2 = 0
      NSURF = 0
      NSUNEW = 0
      DO I=1,IVOL
        IEE = 1+(I-1)/MODTET
        IF (ITOUS.EQ.0.OR.ISDVU(ISDEL(IEE)).GT.0) THEN
          IBB = IELBOR(IEE)
          MODIM1 = 1+MOD(I-1,5)
          IG5(MODIM1) = 0
          VVMA = MAX(VALX(NTET(1,I)),VALX(NTET(2,I))
     &              ,VALX(NTET(3,I)),VALX(NTET(4,I)))
          IF (VVMA.GE.VISOMI) THEN
            VVMI = MIN(VALX(NTET(1,I)),VALX(NTET(2,I))
     &                ,VALX(NTET(3,I)),VALX(NTET(4,I)))
            IF (VVMI.LE.VISOMA) THEN
              IG5(MODIM1) = 1
              NSURF = NSURF+1
              IF (NSURF.GT.NFMAX) THEN
                CALL ERRSURF(ILANG,NFMAX,NSURF,0)
                RETURN
              ENDIF
              NFIS(5,NSURF) = 0
              NEIS(NSURF) = IEE
c
c effet de bord
c
              NBB = 0
              K = 0
              DO WHILE(K.LT.4.AND.NBB.EQ.0)
                K = K+1
                I1 = NTET(IFA(1,K),I)
                I2 = NTET(IFA(2,K),I)
                I3 = NTET(IFA(3,K),I)
                VVMA = MAX(VALX(I1),VALX(I2),VALX(I3))
                VVMI = MIN(VALX(I1),VALX(I2),VALX(I3))
                IF (VVMI.GT.VISOMI.AND.VVMA.LT.VISOMA) THEN
                  NBB = K
                  KMI = KCOMPL(K)
                  XIS(1,NSURF) = X(I1)
                  XIS(2,NSURF) = X(I2)
                  XIS(3,NSURF) = X(I3)
                  YIS(1,NSURF) = Y(I1)
                  YIS(2,NSURF) = Y(I2)
                  YIS(3,NSURF) = Y(I3)
                  ZIS(1,NSURF) = Z(I1)
                  ZIS(2,NSURF) = Z(I2)
                  ZIS(3,NSURF) = Z(I3)
                  XIS(4,NSURF) = XIS(3,NSURF)
                  YIS(4,NSURF) = YIS(3,NSURF)
                  ZIS(4,NSURF) = ZIS(3,NSURF)
                  IF (IFVISO.NE.0) THEN
                    VIS(1,NSURF) = VALXB(I1)
                    VIS(2,NSURF) = VALXB(I2)
                    VIS(3,NSURF) = VALXB(I3)
                    VIS(4,NSURF) = VIS(3,NSURF)
                  ENDIF
               ENDIF
              ENDDO
              NBMI = 0
              NTP = 0
              ITRIV = 0
              IF (NBB.EQ.0) THEN
                DO K=1,4
                  IF (VALX(NTET(K,I)).LE.VISO) THEN
                    NBMI = NBMI+1
                    KMI = NTET(K,I)
                  ENDIF
                ENDDO
                IF (NBMI.GT.0.AND.NBMI.LT.4) THEN
                  IF (NBMI.EQ.2) THEN
                    NINTER = 4
                  ELSE
                    NINTER = 3
                  ENDIF
                  DO K=1,6
                    IF (NINTER.GT.0) THEN
                      N1 = NTET(IAR(1,K),I)
                      N2 = NTET(IAR(2,K),I)
                      V1 = MAX(VALX(N1),VALX(N2))
                      V2 = MIN(VALX(N1),VALX(N2))
                      IF (VISO.LE.V1.AND.VISO.GE.V2) THEN
                        IF (V1.EQ.V2) THEN
                          IF (IFVISO.NE.0) THEN
                            VIS(NINTER,NSURF)   = VALXB(N1)
                            VIS(NINTER-1,NSURF) = VALXB(N2)
                          ENDIF
                          XIS(NINTER,NSURF) = X(N1)
                          YIS(NINTER,NSURF) = Y(N1)
                          ZIS(NINTER,NSURF) = Z(N1)
                          NINTER = NINTER-1
                          IF (NINTER.GT.0) THEN
                            XIS(NINTER,NSURF) = X(N2)
                            YIS(NINTER,NSURF) = Y(N2)
                            ZIS(NINTER,NSURF) = Z(N2)
                            NINTER = NINTER-1
                          ENDIF
                        ELSE
                          T = (V1-VISO)/(V1-V2)
                          IF (T.LT.EPST.OR.T.GT.UMEPST) NTP = NTP+1
                          IF (V1.NE.VALX(N1)) THEN
                            NBID = N1
                            N1 = N2
                            N2 = NBID
                          ENDIF 
                          IF (IFVISO.NE.0)
     &                         VIS(NINTER,NSURF) = T*VALXB(N2)
     &                                      + (1.-T)*VALXB(N1)
                          XIS(NINTER,NSURF) = T*X(N2) + (1.-T)*X(N1)
                          YIS(NINTER,NSURF) = T*Y(N2) + (1.-T)*Y(N1)
                          ZIS(NINTER,NSURF) = T*Z(N2) + (1.-T)*Z(N1)
                          NINTER = NINTER-1
                        ENDIF
                      ENDIF
                    ENDIF
                  ENDDO
ccc            IF (NINTER.NE.0) PRINT*,'bizarre et etrange',NSURF
                 IF (NBMI.EQ.2) THEN
                    CALL CONVEX(XIS(1,NSURF),YIS(1,NSURF),ZIS(1,NSURF)
     &                         ,VIS(1,NSURF),IFVISO)
                  ELSE
                    XIS(4,NSURF) = XIS(3,NSURF)
                    YIS(4,NSURF) = YIS(3,NSURF)
                    ZIS(4,NSURF) = ZIS(3,NSURF)
                    IF (IFVISO.NE.0) VIS(4,NSURF) = VIS(3,NSURF)
                  ENDIF
                  CALL ELIMINDEGEN(XIS(1,NSURF),YIS(1,NSURF)
     &                            ,ZIS(1,NSURF),ITRIV)
                  IF (NTP.GT.0) ITRIV = 1
c              IF (ITRIV.NE.0) print*,'*** elimine'
                ELSE
                  ITRIV = 1
                ENDIF
              ENDIF
C
              IF (ITRIV.NE.0) THEN
                NSURF = NSURF-1
              ELSE
C
C On met les faces dans le bon sens
C
                IF (ICALSU.GE.0.AND.NBB.EQ.0) THEN
                  PMIXT = (XIS(2,NSURF)-XIS(1,NSURF))*(
     &               (YIS(3,NSURF)-YIS(1,NSURF))*(Z(KMI)-ZIS(1,NSURF))
     &              -(ZIS(3,NSURF)-ZIS(1,NSURF))*(Y(KMI)-YIS(1,NSURF)))
     &                - (YIS(2,NSURF)-YIS(1,NSURF))*(
     &               (XIS(3,NSURF)-XIS(1,NSURF))*(Z(KMI)-ZIS(1,NSURF))
     &              -(ZIS(3,NSURF)-ZIS(1,NSURF))*(X(KMI)-XIS(1,NSURF)))
     &                + (ZIS(2,NSURF)-ZIS(1,NSURF))*(
     &               (XIS(3,NSURF)-XIS(1,NSURF))*(Y(KMI)-YIS(1,NSURF))
     &              -(YIS(3,NSURF)-YIS(1,NSURF))*(X(KMI)-XIS(1,NSURF)))
                  IF ((PMIXT.GT.0..AND.ICALSU.LE.1)
     &            .OR.(PMIXT.LT.0..AND.ICALSU.EQ.2)) THEN
                    IF (NBMI.EQ.2) THEN
                      CALL ECHR(XIS(2,NSURF),XIS(4,NSURF))
                      CALL ECHR(YIS(2,NSURF),YIS(4,NSURF))
                      CALL ECHR(ZIS(2,NSURF),ZIS(4,NSURF))
                      IF (IFVISO.NE.0)
     &                     CALL ECHR(VIS(2,NSURF),VIS(4,NSURF))
                    ELSE
                      CALL ECHR(XIS(2,NSURF),XIS(3,NSURF))
                      CALL ECHR(YIS(2,NSURF),YIS(3,NSURF))
                      CALL ECHR(ZIS(2,NSURF),ZIS(3,NSURF))
                      XIS(4,NSURF) = XIS(3,NSURF)
                      YIS(4,NSURF) = YIS(3,NSURF)
                      ZIS(4,NSURF) = ZIS(3,NSURF)
                      IF (IFVISO.NE.0) THEN
                        CALL ECHR(VIS(2,NSURF),VIS(3,NSURF))
                        VIS(4,NSURF) = VIS(3,NSURF)
                      ENDIF
                    ENDIF
                  ENDIF
                ENDIF
                IF (MOD(I,5).NE.0.AND.IBB.NE.0.AND.ICALSU.GT.0) THEN
                  NSURF0 = NSURF
                  DO K=2,4
                    I1 = NTET(IFA(1,K),I)
                    I2 = NTET(IFA(2,K),I)
                    I3 = NTET(IFA(3,K),I)
                    VVMA = MAX(VALX(I1),VALX(I2),VALX(I3))
                    VVMI = MIN(VALX(I1),VALX(I2),VALX(I3))
                    IF (VVMA.GT.VISOmi.AND.VVMI.LT.VISOma) THEN
                      IF (IBB.LT.0) THEN
                        ICOM = 0
                        DO M=1,NF
                          ICO = 0
                          DO L=1,NDS
                            IF (I1.EQ.NFAC(L,M).OR.I2.EQ.NFAC(L,M)
     &                      .OR.I3.EQ.NFAC(L,M)) ICO = ICO+1
                          ENDDO
                          ICOM = MAX(ICOM,ICO)
                          IF (ICOM.EQ.3) GOTO 4444
                        ENDDO
                      ELSE
                        ICO = 0
                        DO L=1,NDS
                          IF (I1.EQ.NFAC(L,IBB).OR.I2.EQ.NFAC(L,IBB)
     &                    .OR.I3.EQ.NFAC(L,IBB)) ICO = ICO+1
                        ENDDO
                        IF (ICO.EQ.3) GOTO 4444
                      ENDIF
                      GOTO 3333
 4444                 IEPS = 0
                      IBUGMI = 10
                      IEPSMI = 0
 5555                 IEPS = IEPS+1
                      EPS = EPSBES(IEPS)
 5556                 EPS2 = EPS**2
                      NBIZ2LO = 0
                      PV1 = (Y(I2)-Y(I3))*(Z(I3)-Z(I1))
     &                    - (Z(I2)-Z(I3))*(Y(I3)-Y(I1))
                      PV2 = (Z(I2)-Z(I3))*(X(I3)-X(I1))
     &                    - (X(I2)-X(I3))*(Z(I3)-Z(I1))
                      PV3 = (X(I2)-X(I3))*(Y(I3)-Y(I1))
     &                    - (Y(I2)-Y(I3))*(X(I3)-X(I1))
                      PVN = PV1**2 + PV2**2 + PV3**2
                      IF (PVN.LT.EPS2) GOTO 3333
                      USPVN = 1./SQRT(PVN)
                      PV1 = PV1*USPVN
                      PV2 = PV2*USPVN
                      PV3 = PV3*USPVN
                      NBINT = 0
                      DO J=1,4
                        JJ = ISUC(J)
                        V1 = XIS(JJ,NSURF0)-XIS(J,NSURF0)
                        V2 = YIS(JJ,NSURF0)-YIS(J,NSURF0)
                        V3 = ZIS(JJ,NSURF0)-ZIS(J,NSURF0)
                        VNOR = V1**2 + V2**2 + V3**2
                        IF (VNOR.GE.EPS2) THEN
                          PMIXT = (PV1*V1 + PV2*V2 + PV3*V3)/SQRT(VNOR)
                          IF (ABS(PMIXT).LT.EPS) THEN
                            V1 = X(I1)-XIS(J,NSURF0)
                            V2 = Y(I1)-YIS(J,NSURF0)
                            V3 = Z(I1)-ZIS(J,NSURF0)
                            VNOR = V1**2 + V2**2 + V3**2
                            IF (VNOR.LT.EPS2) THEN
                              V1 = X(I2)-XIS(J,NSURF0)
                              V2 = Y(I2)-YIS(J,NSURF0)
                              V3 = Z(I2)-ZIS(J,NSURF0)
                              VNOR = V1**2 + V2**2 + V3**2
                              IF (VNOR.LT.EPS2) THEN
                                V1 = X(I3)-XIS(J,NSURF0)
                                V2 = Y(I3)-YIS(J,NSURF0)
                                V3 = Z(I3)-ZIS(J,NSURF0)
                                VNOR = V1**2 + V2**2 + V3**2
                              ENDIF
                            ENDIF
                            PMIXT = MAX(PMIXT
     &                           ,(PV1*V1 + PV2*V2 + PV3*V3)/SQRT(VNOR))
                            IF (ABS(PMIXT).LT.EPS) THEN
                              IF (NBINT.EQ.1) NBIZ2LO = NBIZ2LO+1
                              NBINT = NBINT+1
                              XTMP(1) = XIS(JJ,NSURF0)
                              YTMP(1) = YIS(JJ,NSURF0)
                              ZTMP(1) = ZIS(JJ,NSURF0)
                              XTMP(2) = XIS(J,NSURF0)
                              YTMP(2) = YIS(J,NSURF0)
                              ZTMP(2) = ZIS(J,NSURF0)
                              IF (IFVISO.NE.0) THEN
                                VTMP(1) = VIS(JJ,NSURF0)
                                VTMP(2) = VIS(J,NSURF0)
                              ENDIF
                            ENDIF
                          ENDIF
                        ENDIF
                      ENDDO
                      NBMI = 0
                      IF (NBINT.NE.1) THEN
                        IF (IEPS.LT.NUMEPS) THEN
                          IF (NBIZ2LO+1.LT.IBUGMI) THEN
                            IBUGMI = NBIZ2LO+1
                            IEPSMI = IEPS
                          ENDIF
                          IEPS = IEPS+1
                          GOTO 5555
                        ELSEIF(IEPS.EQ.NUMEPS) THEN
                          EPS = EPSBES(IEPSMI)
                          IEPS = NUMEPS+1
                          GOTO 5556
                        ELSE
                          NBIZ0 = NBIZ0+1
                          NBIZ2 = NBIZ2+NBIZ2LO
                          GOTO 3333
                        ENDIF
                      ENDIF
                      IF (ICALSU.EQ.1) THEN
                        DO J=1,3
                          IF (VALX(NTET(IFA(J,K),I)).LE.VISOma) THEN
                            NBMI = NBMI+1
                            ITMP(NBMI) = NTET(IFA(J,K),I)
                          ENDIF
                        ENDDO
                      ELSE
                        DO J=1,3
                          IF (VALX(NTET(IFA(J,K),I)).GE.VISOmi) THEN
                            NBMI = NBMI+1
                            ITMP(NBMI) = NTET(IFA(J,K),I)
                          ENDIF
                        ENDDO
                      ENDIF
                      NSURF = NSURF+1
                      IF (NSURF.GT.NFMAX) THEN
                        CALL ERRSURF(ILANG,NFMAX,NSURF,0)
                        RETURN
                      ENDIF
                      NEIS(NSURF) = IEE
                      NFIS(5,NSURF) = ICALSU
                      XIS(1,NSURF) = XTMP(1)
                      YIS(1,NSURF) = YTMP(1)
                      ZIS(1,NSURF) = ZTMP(1)
                      XIS(2,NSURF) = XTMP(2)
                      YIS(2,NSURF) = YTMP(2)
                      ZIS(2,NSURF) = ZTMP(2)
                      XIS(3,NSURF) = X(ITMP(1))
                      YIS(3,NSURF) = Y(ITMP(1))
                      ZIS(3,NSURF) = Z(ITMP(1))
                      IF (IFVISO.NE.0) THEN
                        VIS(1,NSURF) = VTMP(1)
                        VIS(2,NSURF) = VTMP(2)
                        VIS(3,NSURF) = VALXB(ITMP(1))
                      ENDIF
                      IF (NBMI.EQ.1) THEN
                        XIS(4,NSURF) = XIS(3,NSURF)
                        YIS(4,NSURF) = YIS(3,NSURF)
                        ZIS(4,NSURF) = ZIS(3,NSURF)
                        IF (IFVISO.NE.0) VIS(4,NSURF) = VIS(3,NSURF)
                      ELSE
                        XIS(4,NSURF) = X(ITMP(2))
                        YIS(4,NSURF) = Y(ITMP(2))
                        ZIS(4,NSURF) = Z(ITMP(2))
                        IF (IFVISO.NE.0) VIS(4,NSURF) = VALXB(ITMP(2))
                        CALL CONVEX(XIS(1,NSURF),YIS(1,NSURF)
     &                             ,ZIS(1,NSURF),VIS(1,NSURF),IFVISO)
                      ENDIF
                      CALL ELIMINDEGEN(XIS(1,NSURF),YIS(1,NSURF)
     &                                ,ZIS(1,NSURF),ITRIV)
                      IF (ITRIV.NE.0) NSURF = NSURF-1
                    ELSEIF(IBB.NE.0.AND.NDS.EQ.4
     &                 .AND.((ICALSU.EQ.1.AND.VVMA.LE.VISOma)
     &                  .OR. (ICALSU.EQ.2.AND.VVMI.GE.VISOmi))) THEN
                      IF (IBB.GT.0) THEN
                        JJ1 = 0
                        JJ2 = 0
                      ELSE
                        JJ1 = -5
                        JJ2 = 5
                      ENDIF
                      DO JJ=JJ1,JJ2
                        ICO = 0
                        IBJ = IABS(IBB)+JJ
                        DO LL=1,4
                          IF (I1.EQ.NFAC(LL,IBJ).OR.I2.EQ.NFAC(LL,IBJ)
     &                    .OR.I3.EQ.NFAC(LL,IBJ)) ICO = ICO+1
                        ENDDO
                        IF (ICO.EQ.3) THEN
                          IF (ICALSU.EQ.1) THEN
                            VVMA = MAX(VALX(NFAC(1,IBJ))
     &                                ,VALX(NFAC(2,IBJ))
     &                                ,VALX(NFAC(3,IBJ))
     &                                ,VALX(NFAC(4,IBJ)))
                            IF (VVMA.GT.VISO) THEN
                              IFLAG = 1
                            ELSE
                              IFLAG = 0
                            ENDIF
                          ELSE
                            VVMI = MIN(VALX(NFAC(1,IBJ))
     &                                ,VALX(NFAC(2,IBJ))
     &                                ,VALX(NFAC(3,IBJ))
     &                                ,VALX(NFAC(4,IBJ)))
                            IF (VVMI.LT.VISO) THEN
                              IFLAG = 1
                            ELSE
                              IFLAG = 0
                            ENDIF
                          ENDIF
                          IF (IFLAG.EQ.1) THEN
                            NSUNEW = NSUNEW+1
                            NSURF = NSURF+1
                            IF (NSURF.GT.NFMAX) THEN
                              CALL ERRSURF(ILANG,NFMAX,NSURF,0)
                              RETURN
                            ENDIF
                            NEIS(NSURF) = IEE
                            NFIS(5,NSURF) = ICALSU
                            IF (IFVISO.NE.0) THEN
                              DO J=1,3
                                XIS(J,NSURF) = X(NTET(IFA(J,K),I))
                                YIS(J,NSURF) = Y(NTET(IFA(J,K),I))
                                ZIS(J,NSURF) = Z(NTET(IFA(J,K),I))
                                VIS(J,NSURF) = VALXB(NTET(IFA(J,K),I))
                              ENDDO
                              XIS(4,NSURF) = XIS(3,NSURF)
                              YIS(4,NSURF) = YIS(3,NSURF)
                              ZIS(4,NSURF) = ZIS(3,NSURF)
                              VIS(4,NSURF) = VIS(3,NSURF)
                            ELSE
                              DO J=1,3
                                XIS(J,NSURF) = X(NTET(IFA(J,K),I))
                                YIS(J,NSURF) = Y(NTET(IFA(J,K),I))
                                ZIS(J,NSURF) = Z(NTET(IFA(J,K),I))
                              ENDDO
                              XIS(4,NSURF) = XIS(3,NSURF)
                              YIS(4,NSURF) = YIS(3,NSURF)
                              ZIS(4,NSURF) = ZIS(3,NSURF)
                            ENDIF
                            CALL ELIMINDEGEN(XIS(1,NSURF),YIS(1,NSURF)
     &                                      ,ZIS(1,NSURF),ITRIV)
                            IF (ITRIV.NE.0) NSURF = NSURF-1
                          ENDIF
                        ENDIF
                      ENDDO
                    ENDIF
 3333               CONTINUE
                  ENDDO
                ENDIF
              ENDIF
            ENDIF
          ENDIF
          IF (MOD(I,5).EQ.0.AND.ICALSU.NE.0.AND.IBB.NE.0) THEN
            IF (MAX(IG5(1),IG5(2),IG5(3),IG5(4)).GT.0) THEN
              DO K=1,4
                IF (IG5(5-K).EQ.0) THEN
                  VV = VALX(NTET(1,I-K))
                  IF ((ICALSU.EQ.1.AND.VV.LE.VISOma)
     &            .OR.(ICALSU.EQ.2.AND.VV.GE.VISOmi)) THEN
                    DO L=1,4
                      I1 = NTET(IFA(1,L),I-K)
                      I2 = NTET(IFA(2,L),I-K)
                      I3 = NTET(IFA(3,L),I-K)
                      IF (IBB.GT.0) THEN
                        ICO = 0
                        DO LL=1,NDS
                          IF (I1.EQ.NFAC(LL,IBB).OR.I2.EQ.NFAC(LL,IBB)
     &                    .OR.I3.EQ.NFAC(LL,IBB)) ICO = ICO+1
                        ENDDO
                      ELSE
                        ICO = 3
                      ENDIF
                      IF (ICO.EQ.3) THEN
                        NSURF = NSURF+1
                        IF (NSURF.GT.NFMAX) THEN
                          CALL ERRSURF(ILANG,NFMAX,NSURF,0)
                          RETURN
                        ENDIF
                        NEIS(NSURF) = IEE
                        NFIS(5,NSURF) = ICALSU
                        IF (IFVISO.NE.0) THEN
                          DO J=1,3
                            XIS(J,NSURF) = X(NTET(IFA(J,L),I-K))
                            YIS(J,NSURF) = Y(NTET(IFA(J,L),I-K))
                            ZIS(J,NSURF) = Z(NTET(IFA(J,L),I-K))
                            VIS(J,NSURF) = VALXB(NTET(IFA(J,L),I-K))
                          ENDDO
                          XIS(4,NSURF) = XIS(3,NSURF)
                          YIS(4,NSURF) = YIS(3,NSURF)
                          ZIS(4,NSURF) = ZIS(3,NSURF)
                          VIS(4,NSURF) = VIS(3,NSURF)
                        ELSE
                          DO J=1,3
                            XIS(J,NSURF) = X(NTET(IFA(J,L),I-K))
                            YIS(J,NSURF) = Y(NTET(IFA(J,L),I-K))
                            ZIS(J,NSURF) = Z(NTET(IFA(J,L),I-K))
                          ENDDO
                          XIS(4,NSURF) = XIS(3,NSURF)
                          YIS(4,NSURF) = YIS(3,NSURF)
                          ZIS(4,NSURF) = ZIS(3,NSURF)
                        ENDIF
                        CALL ELIMINDEGEN(XIS(1,NSURF),YIS(1,NSURF)
     &                                  ,ZIS(1,NSURF),ITRIV)
                        IF (ITRIV.NE.0) NSURF = NSURF-1
                      ENDIF
                    ENDDO
                  ENDIF
                ENDIF
              ENDDO
            ENDIF
          ENDIF
        ENDIF
      ENDDO
C
Cfj      IF (IDERN.EQ.0.AND.ICALSU.NE.0) THEN
Cfj        IF (IPASS.LT.NUMEPS.AND.(NBIZ0+NBIZ2).NE.0) THEN
Cfj          IF ((NBIZ0+NBIZ2).LE.NBIZB) THEN
Cfj            IBEST = IPASS
Cfj            NBIZB = NBIZ0+NBIZ2
Cfj          ENDIF
Cfj          GOTO 1
Cfj        ELSEIF(IPASS.EQ.NUMEPS) THEN
Cfj          IDERN = 1
Cfj          IF ((NBIZ0+NBIZ2).LE.NBIZB) THEN
Cfj            IBEST=IPASS
Cfj          ELSE
Cfj            EPS = EPSBES(IBEST)
Cfj          ENDIF
Cfj          PRINT*,'Meilleur score avec eps =',EPS
Cfj          EPSISO = EPS
Cfj          IF (IBEST.NE.IPASS) GOTO 2
Cfj        ELSE
Cfj          EPSISO = EPS
Cfj          PRINT*,'eps ok =',EPS
Cfj        ENDIF
Cfj      ENDIF
C
      IF (ISTDOUT.EQ.0.AND.IMP.GT.0) THEN
        IF (NBIZ0.GT.0) PRINT*,'*** calsur',NBIZ0,' bug 0'
        IF (NBIZ2.GT.0) PRINT*,'*** calsur',NBIZ2,' bug 2'
C
        IF (ILANG.EQ.0) THEN
          PRINT*,'nsurf =',NSURF,' pour',IVOL,' tétraèdres'
        ELSE
          PRINT*,'nsurf =',NSURF,' for',IVOL,' tetraedrons'
        ENDIF
      ENDIF
C
      NNNN0 = NSURF
      IF (ICALSU.EQ.1) THEN
        DO I=1,NF
          NSD = MOD(ISD2(I),1000)
          IF (ITOUS.EQ.0.OR.ISDVU(NSD).GT.0) THEN
            IF (NDS.EQ.3) THEN
              VVMA = MAX(VALX(NFAC(1,I)),VALX(NFAC(2,I))
     &                  ,VALX(NFAC(3,I)))
            ELSE
              VVMA = MAX(VALX(NFAC(1,I)),VALX(NFAC(2,I))
     &                  ,VALX(NFAC(3,I)),VALX(NFAC(4,I)))
            ENDIF
            IF (VVMA.LE.VISOma) THEN
              NSURF = NSURF+1
              IF (NSURF.GT.NFMAX) THEN
                CALL ERRSURF(ILANG,NFMAX,NSURF,1)
                GOTO 1111
              ENDIF
              NEIS(NSURF) = NNUMFA(I)
              NFIS(5,NSURF) = ICALSU
              IF (IFVISO.NE.0) THEN
                DO K=1,NDS
                  XIS(K,NSURF) = X(NFAC(K,I))
                  YIS(K,NSURF) = Y(NFAC(K,I))
                  ZIS(K,NSURF) = Z(NFAC(K,I))
                  VIS(K,NSURF) = VALXB(NFAC(K,I))
                ENDDO
                IF (NDS.EQ.3) THEN
                  XIS(4,NSURF) = XIS(3,NSURF)
                  YIS(4,NSURF) = YIS(3,NSURF)
                  ZIS(4,NSURF) = ZIS(3,NSURF)
                  VIS(4,NSURF) = VIS(3,NSURF)
                ENDIF
              ELSE
                DO K=1,NDS
                  XIS(K,NSURF) = X(NFAC(K,I))
                  YIS(K,NSURF) = Y(NFAC(K,I))
                  ZIS(K,NSURF) = Z(NFAC(K,I))
                ENDDO
                IF (NDS.EQ.3) THEN
                  XIS(4,NSURF) = XIS(3,NSURF)
                  YIS(4,NSURF) = YIS(3,NSURF)
                  ZIS(4,NSURF) = ZIS(3,NSURF)
                ENDIF
              ENDIF
            ENDIF
          ENDIF
        ENDDO
 1111   IF (ISTDOUT.EQ.0.AND.IMP.GT.0) THEN
          IF (ILANG.EQ.0) THEN
            PRINT*,'nsurf =',NSURF
     &           ,'  (',NSURF-NNNN0,' rajoutées 2ième phase)'
          ELSE
            PRINT*,'nsurf =',NSURF
     &           ,'  (',NSURF-NNNN0,' added 2nd phase)'
          ENDIF
        ENDIF
      ELSEIF(ICALSU.EQ.2) THEN
        DO I=1,NF
          NSD = MOD(ISD2(I),1000)
          IF (ITOUS.EQ.0.OR.ISDVU(NSD).GT.0) THEN
            IF (NDS.EQ.3) THEN
              VVMI = MIN(VALX(NFAC(1,I)),VALX(NFAC(2,I))
     &                  ,VALX(NFAC(3,I)))
            ELSE
              VVMI = MIN(VALX(NFAC(1,I)),VALX(NFAC(2,I))
     &                  ,VALX(NFAC(3,I)),VALX(NFAC(4,I)))
            ENDIF
            IF (VVMI.GE.VISOmi) THEN
              NSURF = NSURF+1
              IF (NSURF.GT.NFMAX) THEN
                CALL ERRSURF(ILANG,NFMAX,NSURF,1)
                GOTO 2222
              ENDIF
              NEIS(NSURF) = NNUMFA(I)
              NFIS(5,NSURF) = ICALSU
              IF (IFVISO.NE.0) THEN
                DO K=1,NDS
                  XIS(K,NSURF) = X(NFAC(K,I))
                  YIS(K,NSURF) = Y(NFAC(K,I))
                  ZIS(K,NSURF) = Z(NFAC(K,I))
                  VIS(K,NSURF) = VALXB(NFAC(K,I))
                ENDDO
                IF (NDS.EQ.3) THEN
                  XIS(4,NSURF) = XIS(3,NSURF)
                  YIS(4,NSURF) = YIS(3,NSURF)
                  ZIS(4,NSURF) = ZIS(3,NSURF)
                  VIS(4,NSURF) = VIS(3,NSURF)
                ENDIF
              ELSE
                DO K=1,NDS
                  XIS(K,NSURF) = X(NFAC(K,I))
                  YIS(K,NSURF) = Y(NFAC(K,I))
                  ZIS(K,NSURF) = Z(NFAC(K,I))
                ENDDO
                IF (NDS.EQ.3) THEN
                  XIS(4,NSURF) = XIS(3,NSURF)
                  YIS(4,NSURF) = YIS(3,NSURF)
                  ZIS(4,NSURF) = ZIS(3,NSURF)
                ENDIF
              ENDIF
            ENDIF
          ENDIF
        ENDDO
 2222   IF (ISTDOUT.EQ.0.AND.IMP.GT.0) THEN
          IF (ILANG.EQ.0) THEN
            PRINT*,'nsurf =',NSURF
     &         ,'  (',NSURF-NNNN0,' rajoutées 2ième phase) (',NSUNEW,')'
          ELSE
            PRINT*,'nsurf =',NSURF
     &           ,'  (',NSURF-NNNN0,' added 2nd phase) (',NSUNEW,')'
          ENDIF
        ENDIF
      ENDIF
      IF (IMP.GT.0) THEN
        CALL TEMPS(S1,IS)
        IF (ISTDOUT.EQ.0.AND.(S1-S0).GE.1.) THEN
          IF (ILANG.EQ.0) THEN
            PRINT*,'temps:',REAL(S1-S0)
          ELSE
            PRINT*,'time:',REAL(S1-S0)
          ENDIF
        ENDIF
      ENDIF
C
      NSURF0 = NSURF
      CALL SYMISO(NSURF0,NSURF,XIS,YIS,ZIS,VIS,NFIS,NEIS
     &           ,NFMAX,NRECON,XMED,YMED,0)
      IF (NSURF.NE.NSURF0) THEN
        IF (ISTDOUT.EQ.0.AND.IMP.GT.0) THEN
          IF (ILANG.EQ.0) THEN
            PRINT*,'nsurf =',NSURF,' après symétries'
          ELSE
            PRINT*,'nsurf =',NSURF,' after symmetries'
          ENDIF
        ENDIF
      ENDIF
C
      CALL INV3X3(ROTA,ROTLOC,IERR)
      DO I=1,3
        DO J=1,3
          ROTA(J,I) = ROTLOC(J,I)
        ENDDO
      ENDDO
      CALL ROTATE(0)
      IELISO = 0
      IF (ICTFAC.GT.15.AND.ICTFAC.LE.97) THEN
        CALL ELISO(IRC)
        IF (IRC.EQ.0) THEN
          IELISO = 1
        ELSE
          ICTFAC = 99
        ENDIF
      ENDIF
C
      END
C=======================================================================
      SUBROUTINE ELIMINDEGEN(XIS,YIS,ZIS,ITRIV)
      DIMENSION XIS(4),YIS(4),ZIS(4)
C
      NEGAL = 0
      DO K=1,3
        DO L=K+1,4
          DDD = (XIS(K)-XIS(L))**2
     &        + (YIS(K)-YIS(L))**2
     &        + (ZIS(K)-ZIS(L))**2
          IF (DDD.EQ.0.) NEGAL = NEGAL+1
        ENDDO
      ENDDO
      IF (NEGAL.GE.2) THEN
        ITRIV = 1
Cfj        print*,negal
Cfj        do i=1,4
Cfj          print*,xis(i),yis(i),zis(i)
Cfj        enddo
      ELSE
        ITRIV = 0
      ENDIF
      END
C=======================================================================
      SUBROUTINE ERRSURF(ILANG,NFMAX,NSURF,IOPT)
      IF (ILANG.EQ.0) THEN
        PRINT*,'*** Isosurface trop compliquée (nb faces max=',NFMAX,')'
        IF (IOPT.EQ.0) THEN
          PRINT*,'*** Abandon...'
        ELSE
          PRINT*,'*** Abandon 2ième phase...'
        ENDIF
      ELSE
        PRINT*,'*** Isosurface too big (nb faces max=',NFMAX,')'
        IF (IOPT.EQ.0) THEN
          PRINT*,'*** Canceled...'
        ELSE
          PRINT*,'*** Canceled 2nd round...'
        ENDIF
      ENDIF
      IF (IOPT.EQ.0) NSURF = 0
      END
