C QUALI6    SOURCE    GOUNAND   26/07/06    21:15:11     12592          
      SUBROUTINE QUALI6(MELEMX,IELDEB,IELFIN,IMET,IMOMET,XDENS,KCMETR
     $     ,NKPVIR,XVTOL,MLREEL,NDQC,jcritt,pcritq,qcritq)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
C***********************************************************************
C NOM         : QUALI6
C DESCRIPTION : Etant donné un maillage volumique simple MELEMX,
C               on construit la qualité de chacun de ses éléments
C               dans un listreel MLREEL.
C               MELEMX est supposé actif.
C               MLREEL est rendu actif.
C
C     Par rapport à quali2, on utilise xvtol pour mettre
C     le volume d'un élément à 0 s'il est petit.
C     Ceci est important car on utilise le signe pour dégrader la
C     qualité d'un élément (-1)
C
C     Par rapport à quali3, MELEME devient un MELEMX, MLREEL est un
C     segment déjà existant et on introduit les éléments de début et de
C     fin IELDEB et IELFIN qui servent pour MELEMX ET MLREEL (qui sont
C     supposés de même dimension cf. trlver.eso.)
*     Pour MLREEL, comme on ne calcule pas la qualité des éléments
*     contenant le noeud virtuel, on a NDQC qui dit le nombre de qualité
*     calculés (IELDEB sert donc pour MELEMX et MLREEL mais IELFIN
*     uniquement pour MELEMX et NDQC uniquement pour MLREEL).
C
*     Par rapport à quali5, on essaie de simplifier et de regrouper les
*     cas avec, sans métrique + un peu de ménage
*
*     En ce qui concerne jcritq, le critère est généralement le mini de
*     : la qualité, la taille, l'inverse de la taille
*     Suivant le chiffre des centaines c de jcritq :
C     * c=0 : on renvoie le critère complet (mini des 3)
C     * c=1 : on renvoie la qualité
C     * c=2 : on renvoie la taille
C     * c=3 : on renvoie l'inverse de la taille
C
C
C
C LANGAGE     : ESOPE
C AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
C               mél : gounand@semt2.smts.cea.fr
C***********************************************************************
C VERSION    : v1, 27/11/2017, version initiale
C HISTORIQUE : v1, 27/11/2017, création
C HISTORIQUE : v2, 10/12/2025, on met des critères de qualité
C     similaires à DEDUADAP
C HISTORIQUE :
C***********************************************************************
-INC CCGEOME
-INC PPARAM
-INC CCOPTIO
-INC CCREEL
-INC SMCOORD
-INC TMATOP1
*-INC SMETRIQ
      POINTEUR KCMETR.METRIQ
*-INC SMELEMX
-INC SMLREEL
      PARAMETER(NMET=6)
      DIMENSION XMET(NMET)
*      DIMENSION XMET2(2,2)
*      DIMENSION XMET3(3,3)
      DIMENSION XMETV(3,3)
      DIMENSION XJAC(3,3)
      DIMENSION XJTMJ(3,3)
      DIMENSION XMJ(3,3)
      REAL*8 DXI(3)
      REAL*8 DARET(6)
      DIMENSION A(3,3),D(3)
* Derivative of the affine barycentric map
      DIMENSION DABM(4,3)
* Simplex node coordinates
      DIMENSION SNCO(3,4)


      DIMENSION IDXSYM(3,3,3)
      LOGICAL LROT
      LOGICAL LQUAL0
*
* Statement functions
*      DISTA(A,B,C,D)=SQRT((A-C)*(A-C)+(B-D)*(B-D))
*      DISTB(A,B,C,D,E,F)=SQRT((A-D)*(A-D)+(B-E)*(B-E)+(C-F)*(C-F))
*      DISTA(A,B)=SQRT(A*A+B*B)
*      DISTB(A,B,C)=SQRT(A*A+B*B+C*C)
*      DETTRI(A11,A12,A21,A22)=A11*A22-A21*A12
      DETTET(A11,A12,A13,A21,A22,A23,A31,A32,A33)=
     &  A11*(A22*A33-A23*A32)
     &       +A12*(A23*A31-A21*A33)
     &       +A13*(A21*A32-A22*A31)
*
      DATA ((XMETV(I,J),I=1,3),J=1,3) /9*0.D0/
      DATA ((XJAC(I,J),I=1,3),J=1,3)  /9*0.D0/
      DATA ((XJTMJ(I,J),I=1,3),J=1,3) /9*0.D0/
      DATA ((XMJ(I,J),I=1,3),J=1,3) /9*0.D0/
      DATA ((A(I,J),I=1,3),J=1,3) /9*0.D0/
      DATA (D(I),I=1,3) /3*0.D0/
      DATA ((DABM(I,J),I=1,4),J=1,3) /12*0.D0/
      DATA ((SNCO(I,J),I=1,3),J=1,4) /12*0.D0/
      DATA (((IDXSYM(I,J,K),I=1,1),J=1,1),K=1,1) /1/
      DATA (((IDXSYM(I,J,K),I=1,2),J=1,2),K=2,2) /1,2,2,3/
      DATA (((IDXSYM(I,J,K),I=1,3),J=1,3),K=3,3) /1,2,4,2,3,5,4,5,6/
*
*
* Executable statements
*
*  Choix du critere de qualite sans metrique
*  0 : Coupez
*  1 : XALIN2
*  2 : XALIN1
*  3 : 2D : D(2)/D(1)
*      3D : D(3)/D(1)
*     write(ioimp,*) 'quali6: NKPVIR=',NKPVIR
      JCRITC=JCRITT/100
      JCRITQ=JCRITT-(JCRITC*100)
      IF (IIMPI.EQ.666) THEN
         WRITE(IOIMP,*) 'JCRITT,JCRITC,JCRITQ=',JCRITT,JCRITC,JCRITQ
         WRITE(IOIMP,*) 'PCRITQ,QCRITQ=',PCRITQ,QCRITQ
      ENDIF
*
      NDQC=0
      IDIMP1=IDIM+1
*      NBNN=NUMX(/1)
      NBNN=NNCOU
*      NBELEM=NUMX(/2)
      IF (NBNN.NE.IDIMP1) THEN
         CALL ERREUR(5)
         RETURN
      ENDIF
      IF
     $     (.NOT.(IELDEB.GE.1.AND.IELFIN.GE.IELDEB.AND.NLCOU.GE.IELFIN
     $     .AND.NUMX(/2).GE.NLCOU))  THEN
         WRITE(IOIMP,*) 'coucou quali6'
         write(ioimp,*) 'IELDEB=',IELDEB
         write(ioimp,*) 'IELFIN=',IELFIN
         write(ioimp,*) 'NLCOU=',NLCOU
         write(ioimp,*) 'NUM2=',NUMX(/2)
         write(ioimp,*) 'MELEMX=',MELEMX
         CALL ECMELX(MELEMX,0)
         CALL ERREUR(5)
         RETURN
      ENDIF
*         XPET=XPETIT*10.D0
      XPET=sqrt(XPETIT)

*     DO 10 IBELEM=1,NUMX(/2)
      DO 10 IBELEM=IELDEB,IELFIN
         LQUAL0=.FALSE.
*         WRITE(IOIMP,*) 'IBELEM=',IBELEM
*     Calcul de la métrique moyenne M soit dans XMETD (scalaire)
*                                     soit dans XMETV (tenseur SPD)
* Derivative of the affine barycentric map M : lambda(x) = Mx +c
*     Les coordonnees barycentriques sont definies par rapport au
*     simplex regulier de cote 1, centre sur l'origine. Le noeud sommet
*     a toutes ses coordonnees nulles sauf la derniere
* Initialisations au premier pas
         IF (IBELEM.EQ.IELDEB) THEN
            IF (IDIM.GE.1) THEN
               DABM(1,1)=-1.D0
               DABM(2,1)=+1.D0
               IF (IDIM.GE.2) THEN
                  DABM(1,2)=-1.D0/SQRT(3.D0)
                  DABM(2,2)=-1.D0/SQRT(3.D0)
                  DABM(3,2)=+2.D0/SQRT(3.D0)
                  IF (IDIM.GE.3) THEN
                     DABM(1,3)=-1.D0/SQRT(6.D0)
                     DABM(2,3)=-1.D0/SQRT(6.D0)
                     DABM(3,3)=-1.D0/SQRT(6.D0)
                     DABM(4,3)=+3.D0/SQRT(6.D0)
                  ENDIF
               ENDIF
            ENDIF
*
            IF (IMET.EQ.1) XMETD=1.D0/DENSIT
            IF (IMET.EQ.2) XMETD=1.D0/XDENS
            IF (IMET.EQ.3) NFMET=1
            IF (IMET.EQ.4) NFMET=IDIM*(IDIM+1)/2
         ENDIF
*
         if (imet.gt.0) then
            IF (IMET.GE.1.AND.IMET.LE.3) THEN
               IF (IMET.EQ.3) THEN
                  YDENS=0.D0
                  DO I=1,IDIMP1
                     INO=NUMX(I,IBELEM)
                     IF (NKPVIR.NE.0) THEN
                        IF (INO.LE.NKPVIR) GOTO 10
                     ENDIF
*     Ici on fait la moyenne arithmétique
*     mais kcmetr contient le log du tenseur si imomet=1
                     YDENS=YDENS+KCMETR.XIN(1,INO)
                  ENDDO
                  YDENS=YDENS/IDIMP1
                  if (imomet.eq.1) then
                     YDENS=EXP(YDENS)
                  endif
                  XMETD2=YDENS
                  XMETD=SQRT(XMETD)
               ELSE
                  XMETD2=XMETD**2
               ENDIF
            ELSEIF (IMET.EQ.4) THEN
               DO J=1,NFMET
                  XMET(J)=0.D0
               ENDDO
               DO I=1,IDIMP1
                  INO=NUMX(I,IBELEM)
                  IF (NKPVIR.NE.0) THEN
                     IF (INO.LE.NKPVIR) GOTO 10
                  ENDIF
                  DO J=1,NFMET
                     XMET(J)=XMET(J)+KCMETR.XIN(J,INO)
                  ENDDO
               ENDDO
               DO J=1,NFMET
                  XMET(J)=XMET(J)/IDIMP1
               ENDDO
*
               if (imomet.eq.1) then
                  DO J=1,IDIM
                     DO I=1,IDIM
                        A(I,J)=XMET(IDXSYM(I,J,IDIM))
                     ENDDO
                  ENDDO
*     Exponentielle du tenseur symétrique
                  IOTENS=8
                  IKAS=3
                  CALL TENS2(IOTENS,IKAS,A,D,XMETV)
                  IF (IERR.NE.0) RETURN
               else
                  DO J=1,IDIM
                     DO I=1,IDIM
                        XMETV(I,J)=XMET(IDXSYM(I,J,IDIM))
                     ENDDO
                  ENDDO
               endif
            ELSE
               WRITE(IOIMP,*) 'quali6 imet=',IMET
               CALL ERREUR(5)
               RETURN
            ENDIF
         endif
* Determinant de la metrique
         if (jcritq.eq.0) then
            if (imet.gt.0) then
               IF (IMET.GE.1.AND.IMET.LE.3) THEN
                  XDETMD=SQRT(XMETD2)**IDIM
               ELSEIF (IMET.EQ.4) THEN
                  IF (IDIM.EQ.1) THEN
                     XDETMD=XMETV(1,1)
                  ELSEIF (IDIM.EQ.2) THEN
                     XDETMD=XMETV(1,1)*XMETV(2,2)-XMETV(2,1)*XMETV(1,2)
                  ELSEIF (IDIM.EQ.3) THEN
                     XDETMD=DETTET(XMETV(1,1),XMETV(1,2),XMETV(1,3)
     $                    ,XMETV(2,1),XMETV(2,2),XMETV(2,3),XMETV(3,1)
     $                    ,XMETV(3,2),XMETV(3,3))
                  ELSE
                     WRITE(IOIMP,*) 'quali6 jcrit=0 imet=4 idim=',IDIM
                     INTERR(1)=IDIM
                     CALL ERREUR(709)
                     RETURN
                  ENDIF
                  XDETMD=SQRT(XDETMD)
               ELSE
                  WRITE(IOIMP,*) 'quali6 imet=',IMET
                  CALL ERREUR(5)
                  RETURN
               ENDIF
            endif
* Volume du simplex
            IF (IDIM.EQ.1) THEN
               I0=NUMX(1,IBELEM)
               I1=NUMX(2,IBELEM)
               IF (NKPVIR.NE.0) THEN
                  IF (I0.LE.NKPVIR.OR.I1.LE.NKPVIR) goto 10
               ENDIF
               IP0=(I0-1)*IDIMP1
               IP1=(I1-1)*IDIMP1
               X10=XCOOR(IP1+1)-XCOOR(IP0+1)
               XVOLO=X10
            ELSEIF (IDIM.EQ.2) THEN
               I0=NUMX(1,IBELEM)
               I1=NUMX(2,IBELEM)
               I2=NUMX(3,IBELEM)
               IF (NKPVIR.NE.0) THEN
                  IF (I0.LE.NKPVIR.OR.I1.LE.NKPVIR.OR.I2.LE.NKPVIR) goto
     $                 10
               ENDIF
               IP0=(I0-1)*IDIMP1
               IP1=(I1-1)*IDIMP1
               IP2=(I2-1)*IDIMP1
               X10=XCOOR(IP1+1)-XCOOR(IP0+1)
               Y10=XCOOR(IP1+2)-XCOOR(IP0+2)
               X20=XCOOR(IP2+1)-XCOOR(IP0+1)
               Y20=XCOOR(IP2+2)-XCOOR(IP0+2)
               XVOLO=(X10*Y20-X20*Y10)/2.D0
            ELSEIF (IDIM.EQ.3) THEN
               I0=NUMX(1,IBELEM)
               I1=NUMX(2,IBELEM)
               I2=NUMX(3,IBELEM)
               I3=NUMX(4,IBELEM)
               IF (NKPVIR.NE.0) THEN
                  IF (I0.LE.NKPVIR.OR.I1.LE.NKPVIR.OR.I2.LE.NKPVIR.OR.
     $                 I3.LE.NKPVIR) goto 10
               ENDIF
               IP0=(I0-1)*IDIMP1
               IP1=(I1-1)*IDIMP1
               IP2=(I2-1)*IDIMP1
               IP3=(I3-1)*IDIMP1
               X10=XCOOR(IP1+1)-XCOOR(IP0+1)
               Y10=XCOOR(IP1+2)-XCOOR(IP0+2)
               Z10=XCOOR(IP1+3)-XCOOR(IP0+3)
               X20=XCOOR(IP2+1)-XCOOR(IP0+1)
               Y20=XCOOR(IP2+2)-XCOOR(IP0+2)
               Z20=XCOOR(IP2+3)-XCOOR(IP0+3)
               X30=XCOOR(IP3+1)-XCOOR(IP0+1)
               Y30=XCOOR(IP3+2)-XCOOR(IP0+2)
               Z30=XCOOR(IP3+3)-XCOOR(IP0+3)
               XVOLO=(DETTET(X10,X20,X30,Y10,Y20,Y30,Z10,Z20,Z30))
     $              /6.D0
            ENDIF
* Déterminant de la métrique
            if (imet.gt.0) then
               XVOLO=XVOLO*XDETMD
*               write(ioimp,*) 'XVOLO,XDETMD=',XVOLO,XDETMD
            endif
            IF (IIMPI.EQ.666) THEN
               write(ioimp,*) 'Xvol,xvolo,xvtol=',Xvol,xvolo,xvtol
            ENDIF
            xvol=abs(xvolo)
            if (xvol.LT.xvtol) then
               LQUAL0=.TRUE.
               goto 666
            endif
* Calcul de la longueur de reference XLAARI
            XLAR=0.D0
            IARET=0
            DO IBNN=1,NBNN
               DO JBNN=IBNN+1,NBNN
                  I0=NUMX(IBNN,IBELEM)
                  I1=NUMX(JBNN,IBELEM)
                  IP0=(I0-1)*IDIMP1
                  IP1=(I1-1)*IDIMP1
                  DO IIDIM=1,IDIM
                     DXI(IIDIM)=XCOOR(IP1+IIDIM)-XCOOR(IP0+IIDIM)
                     IF (IMET.GE.1.AND.IMET.LE.3) THEN
                        DXI(IIDIM)=DXI(IIDIM)*XMETD
                     ENDIF
                  ENDDO
                  DXLAR2=0.D0
                  IF (IMET.EQ.4) THEN
                     DO J=1,IDIM
                        DO I=1,IDIM
                           DXLAR2=DXLAR2+XMETV(I,J)*DXI(I)*DXI(J)
                        ENDDO
                     ENDDO
                  ELSE
                     DO I=1,IDIM
                        DXLAR2=DXLAR2+DXI(I)*DXI(I)
                     ENDDO
                  ENDIF
                  DXLAR=SQRT(DXLAR2)
                  IARET=IARET+1
                  DARET(IARET)=DXLAR
                  XLAR=XLAR+DXLAR
               ENDDO
            ENDDO
            NARET=((NBNN-1)*NBNN)/2
            XLAARI=XLAR/NARET
* Coefficient de normalisation
            I=IDIM
            XCOQ=DFACT(I)*(SQRT((DBLE(2**I))/(DBLE(I+1))))
            XQUALC=((XVOL*XCOQ)**(1.D0/IDIM))/XLAARI
         else
* Calcul du jacobien de la transformation geometrique entre
* l'element regulier de coté 1 et l'element courant
*   Coordonnees des noeuds
            DO J=1,IDIMP1
               INOD=NUMX(J,IBELEM)
               IF (NKPVIR.NE.0) THEN
                  IF (INOD.LE.NKPVIR) goto 10
               ENDIF
               IPNOD=(INOD-1)*IDIMP1
               DO I=1,IDIM
                  SNCO(I,J)=XCOOR(IPNOD+I)
               ENDDO
            ENDDO
*         write(ioimp,*) 'SNCO,I,J=',IDIM,IDIMP1
*         write(ioimp,*) ((SNCO(I,J),I=1,IDIM),J=1,IDIMP1)
*         write(ioimp,*) 'DABM,I,J=',IDIMP1,IDIM
*         write(ioimp,*) ((DABM(I,J),I=1,IDIMP1),J=1,IDIM)
*   Matrice Jacobienne de la transformation J = SNCO*DABM
            DO J=1,IDIM
               DO I=1,IDIM
                  XIJ=0.D0
                  DO K=1,IDIMP1
                     XIJ=XIJ+SNCO(I,K)*DABM(K,J)
                  ENDDO
                  XJAC(I,J)=XIJ
               ENDDO
            ENDDO
            IF (IDIM.EQ.1) THEN
               XDETJ=XJAC(1,1)
            ELSEIF (IDIM.EQ.2) THEN
               XDETJ=XJAC(1,1)*XJAC(2,2)-XJAC(2,1)*XJAC(1,2)
            ELSEIF (IDIM.EQ.3) THEN
               XDETJ=DETTET(XJAC(1,1),XJAC(1,2),XJAC(1,3),XJAC(2
     $              ,1),XJAC(2,2),XJAC(2,3),XJAC(3,1),XJAC(3,2)
     $             ,XJAC(3,3))
            ELSE
               INTERR(1)=IDIM
               CALL ERREUR(709)
               RETURN
            ENDIF
*            write(ioimp,*) 'XDETJ=',XDETJ
*         write(ioimp,*) 'XJAC,I,J=',IDIM,IDIM
*         write(ioimp,*) ((XJAC(I,J),I=1,IDIM),J=1,IDIM)
* Matrice JtMJ
            IF (IMET.LT.4) THEN
               DO K=1,IDIM
                  DO I=1,IDIM
                     XIK=0.D0
                     DO J=1,IDIM
                        XIK=XIK+XJAC(J,I)*XJAC(J,K)
                     ENDDO
                     IF (IMET.EQ.0) THEN
                        XJTMJ(I,K)=XIK
                     ELSE
                        XJTMJ(I,K)=XIK*XMETD2
                     ENDIF
                  ENDDO
               ENDDO
               IF (IMET.EQ.0) THEN
                  XDETM=1.D0
               ELSE
                  XDETM=XMETD2**IDIM
               ENDIF
            ELSE
               DO L=1,IDIM
                  DO J=1,IDIM
                     XJL=0.D0
                     DO K=1,IDIM
* Utilisons la symetrie de XMETV
                        XJL=XJL+XMETV(K,J)*XJAC(K,L)
                     ENDDO
                     XMJ(J,L)=XJL
                  ENDDO
               ENDDO
               DO L=1,IDIM
                  DO I=1,IDIM
                     XIL=0.D0
                     DO J=1,IDIM
                        XIL=XIL+XJAC(J,I)*XMJ(J,L)
                     ENDDO
                     XJTMJ(I,L)=XIL
                  ENDDO
               ENDDO
               IF (IDIM.EQ.1) THEN
                  XDETM=XMETV(1,1)
               ELSEIF (IDIM.EQ.2) THEN
                  XDETM=XMETV(1,1)*XMETV(2,2)-XMETV(2,1)*XMETV(1,2)
               ELSEIF (IDIM.EQ.3) THEN
                  XDETM=DETTET(XMETV(1,1),XMETV(1,2),XMETV(1,3),XMETV(2
     $                 ,1),XMETV(2,2),XMETV(2,3),XMETV(3,1),XMETV(3,2)
     $                 ,XMETV(3,3))
               ELSE
                  INTERR(1)=IDIM
                  CALL ERREUR(709)
                  RETURN
               ENDIF
            ENDIF
*      WRITE(IOIMP,*) 'XDETM=',XDETM
*         write(ioimp,*) 'XJTMJ,I,J=',IDIM,IDIM
*     write(ioimp,*) ((XJTMJ(I,J),I=1,IDIM),J=1,IDIM)
            XDETJM=XDETJ*SQRT(XDETM)
*     WRITE(IOIMP,*) 'XDETJM=',XDETJM
            IF (IIMPI.EQ.666) THEN
               write(ioimp,*) 'XDETJ,XDETM,Xdetjm,xvtol=',xdetj,xdetm
     $              ,Xdetjm,xvtol
            ENDIF
*            if (XDETJM.LT.xvtol) then
            if (ABS(XDETJM).LT.xvtol) then
               LQUAL0=.TRUE.
               goto 666
            endif
* Determinant et trace de JTMJ
            IF (IDIM.EQ.1) THEN
               D(1)=XDETJM**2
            ELSEIF (IDIM.EQ.2) THEN
*     XDET=XJTMJ(1,1)*XJTMJ(2,2)-XJTMJ(2,1)*XJTMJ(1,2)
               CALL JACOD2(XJTMJ,D)
            ELSEIF (IDIM.EQ.3) THEN
*     XDET=DETTET(XJTMJ(1,1),XJTMJ(1,2),XJTMJ(1,3),XJTMJ(2
*     $           ,1),XJTMJ(2,2),XJTMJ(2,3),XJTMJ(3,1),XJTMJ(3,2)
*     $           ,XJTMJ(3,3))
               CALL JACOD3(XJTMJ,3,D)
            ELSE
               WRITE(IOIMP,*) 'quali6 idim=',IDIM
               INTERR(1)=IDIM
               CALL ERREUR(709)
               RETURN
            ENDIF
* On stocke les racines des valeurs propres (longueurs) dans DARET
*         WRITE(IOIMP,*) '1',(D(II),II=1,IDIM)
            DO I=1,IDIM
*     D(I)=ABS(D(I))
               DARET(I)=SQRT(MAX(D(I),XZERO))
            ENDDO
            NARET=IDIM
*
* Pour le calcul de XALIN2 et XALIN1, on stocke D(I) ou DARET(I) dans DXI(I)
*
* Calcul de XALIN2 = inverse de QALI dans Huang
*
*     Calcul de XALIN1 pareil que XALIN2 mais avec les valeurs propres
*     au lieu de leur carré
            IF (JCRITQ.EQ.1.OR.JCRITQ.EQ.2) THEN
               IF (IDIM.EQ.1) THEN
                  XALIN=1
               ELSE
                  if (jcritq.eq.1) then
                     DO I=1,IDIM
                        DXI(I)=D(I)
                     ENDDO
                  else
                     DO I=1,IDIM
                        DXI(I)=DARET(I)
                     ENDDO
                  endif
* Trace
                  XTR=DXI(1)
                  DO I=2,IDIM
                     XTR=XTR+DXI(I)
                  ENDDO
                  XLT=XTR/IDIM
                  if (jcritq.eq.1) then
                     IF (IDIM.EQ.2) THEN
                        XLTD=XLT
                     ELSEIF (IDIM.EQ.3) THEN
                        XLTD=XLT*SQRT(XLT)
                     ELSE
                        INTERR(1)=IDIM
                        CALL ERREUR(709)
                        RETURN
                     ENDIF
                  else
                     XLTD=XLT**IDIM
                  endif
                  XALIN=ABS(XDETJM)/(MAX(XLTD,XPET))
                  IF (IDIM.EQ.3) XALIN=SQRT(XALIN)
               ENDIF
            ENDIF
         endif
* Par rapport au livre de Huang p.205, XALIN2 vaut 1 / (Qali^(n-1)) (n>=2)
* Comme Qali est minore par le rapport d'aspect d'un element, il
* faut comparer XALIN2 a des (rapports de longueur)^(n-1)
* Il y a donc un carre en dimension 3 cf. les expressions de XQUALN
* plus bas (voir aussi deadutil.procedur qui exprime des indicateurs
* en rapports de longueur)
* On a choisi d'exprimer XALIN2 en fonction de XDETJM directement
* car la presque nullite de XDETJ permet de detecter les elements
* plats.
* Si on l'eleve a la puissance (1/3), ca ne marche pas.
* On devra sans doute faire mieux pour etre vraiment robuste....
*         write(ioimp,*) 'XALIN2=',XALIN2
* Les valeurs propres sont censees etre positives mais pas garanti
*     donc on prend la valeur absolue et on reclasse
* Raccourci pour les elements de qualite nulle
 666     CONTINUE
         JELDEB=IELDEB+NDQC
*     write(ioimp,*) 'jeldeb,prog=',jeldeb,prog(jeldeb)
         IF (LQUAL0) THEN
            IF (IIMPI.EQ.666) write(ioimp,*) 'Cas LQUAL0 IELEM=',JELDEB
            PROG(JELDEB)=0.D0
         ELSE
            IF (IDIM.EQ.1) THEN
               XQUALN=1.D0
            ELSE
               IF (JCRITQ.EQ.0) THEN
                  XQUALN=XQUALC
               ELSEIF (JCRITQ.EQ.1.OR.JCRITQ.EQ.2) THEN
                  XQUALN=XALIN
               ELSEIF (JCRITQ.EQ.3) THEN
                  IF (DARET(1).NE.XZERO) THEN
                     XQUALN=DARET(IDIM)/DARET(1)
                  ELSE
                     XQUALN=XZERO
                  ENDIF
               ELSE
                  WRITE(IOIMP,*) 'quali6 : jcritq=',JCRITQ
                  CALL ERREUR(5)
                  RETURN
               ENDIF
            ENDIF
            IF (IMET.EQ.0) THEN
               PROG(JELDEB)=XQUALN
            ELSE
               XQUALQ=XQUALN**(1.d0/QCRITQ)
*     WRITE(IOIMP,*) '3',(DARET(II),II=1,NARET)
               IF (JCRITQ.EQ.0) THEN
                  DMIN=XGRAND
                  DMAX=XZERO
                  DO K=1,NARET
                     DMIN=MIN(DMIN,DARET(K))
                     DMAX=MAX(DMAX,DARET(K))
                  ENDDO
               ELSEIF (JCRITQ.GE.1.AND.JCRITQ.LE.3) THEN
                  DMAX=DARET(1)
                  DMIN=DARET(NARET)
               ELSE
                  WRITE(IOIMP,*) 'quali6 : jcritq=',JCRITQ
                  CALL ERREUR(5)
                  RETURN
               ENDIF
               dinf=abs(DMAX)
               if (dinf.lt.xpet) then
                  XL1=XZERO
               else
                  if (pcritq.eq.0.d0) then
                     DMOYP=1.D0
                     DO K=1,NARET
                        DMOYP=DMOYP*DARET(K)
                     ENDDO
                     XP=1.D0/NARET
                     DMOYP=DMOYP**XP
                  elseif (pcritq.gt.100.d0) then
                     dmoyp=abs(dmax)
                  elseif (pcritq.lt.-100.d0) then
                     dmoyp=abs(dmin)
                  else
                     xsumn=(abs(daret(NARET))/dinf)**pcritq
                     DO K=NARET-1,1,-1
                        xsumn=xsumn+(abs(daret(k))/dinf)**pcritq
                     ENDDO
                     xsumn=xsumn/NARET
                     dmoyp=dinf*xsumn**(1.d0/pcritq)
                  endif
               endif
               XL1=DMOYP
               XL1I=1.D0/MAX(XL1,XPET)
               XL1T=MIN(XL1,XL1I)
*               WRITE(IOIMP,'(A,2X,I3,2X,10(2(f9.3),1X))')
*     $          'critq,d1,d2,dmoyp,xl1t,xalin1,xqual=',
*     $           jcritq,pcritq,qcritq,daret(1),daret(2),dmoyp,xl1t
*     ,xalin1,xqual
               IF (JCRITC.EQ.0) THEN
                  PROG(JELDEB)=MIN(XQUALQ,XL1T)
               ELSEIF (JCRITC.EQ.1) THEN
                  PROG(JELDEB)=XQUALQ
               ELSEIF (JCRITC.EQ.2) THEN
                  PROG(JELDEB)=XL1
               ELSEIF (JCRITC.EQ.3) THEN
                  PROG(JELDEB)=XL1I
               ELSE
                  CALL ERREUR(5)
                  RETURN
               ENDIF
            ENDIF
         ENDIF

* Scaling ? mmmmm, petit doute
         NDQC=NDQC+1
*         Write(ioimp,*) 'jeldeb,prog=',jeldeb,prog(jeldeb)
 10   CONTINUE
      RETURN
*
* formats
*
  188  FORMAT (2X,12(A6,'=',1PG12.5,2X))
*
* End of subroutine QUALI6
*
      END