Télécharger la plaquette Cast3M


fr en

Télécharger quali6.eso

Retour à la liste

Numérotation des lignes : 

quali6

  1. C QUALI6    SOURCE    GOUNAND   26/01/09    21:15:51     12442          
  2.       SUBROUTINE QUALI6(MELEMX,IELDEB,IELFIN,IMET,IMOMET,XDENS,KCMETR
  3.      $     ,NKPVIR,XVTOL,MLREEL,NDQC,ISTRID)
  4.       IMPLICIT REAL*8 (A-H,O-Z)
  5.       IMPLICIT INTEGER (I-N)
  6. C***********************************************************************
  7. C NOM         : QUALI6
  8. C DESCRIPTION : Etant donné un maillage volumique simple MELEMX,
  9. C               on construit la qualité de chacun de ses éléments
  10. C               dans un listreel MLREEL.
  11. C               MELEMX est supposé actif.
  12. C               MLREEL est rendu actif.
  13. C
  14. C     Par rapport à quali2, on utilise xvtol pour mettre
  15. C     le volume d'un élément à 0 s'il est petit.
  16. C     Ceci est important car on utilise le signe pour dégrader la
  17. C     qualité d'un élément (-1)
  18. C
  19. C     Par rapport à quali3, MELEME devient un MELEMX, MLREEL est un
  20. C     segment déjà existant et on introduit les éléments de début et de
  21. C     fin IELDEB et IELFIN qui servent pour MELEMX ET MLREEL (qui sont
  22. C     supposés de même dimension cf. trlver.eso.)
  23. *     Pour MLREEL, comme on ne calcule pas la qualité des éléments
  24. *     contenant le noeud virtuel, on a NDQC qui dit le nombre de qualité
  25. *     calculés (IELDEB sert donc pour MELEMX et MLREEL mais IELFIN
  26. *     uniquement pour MELEMX et NDQC uniquement pour MLREEL).
  27. C
  28. *     Par rapport à quali5, on essaie de simplifier et de regrouper les
  29. *     cas avec, sans métrique + un peu de ménage
  30. *
  31. *     On introduit ISTRID, le nombre d'indicateurs renvoyés
  32. C
  33. C
  34. C
  35. C LANGAGE     : ESOPE
  36. C AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
  37. C               mél : gounand@semt2.smts.cea.fr
  38. C***********************************************************************
  39. C VERSION    : v1, 27/11/2017, version initiale
  40. C HISTORIQUE : v1, 27/11/2017, création
  41. C HISTORIQUE : v2, 10/12/2025, on met des critères de qualité
  42. C     similaires à DEDUADAP
  43. C HISTORIQUE :
  44. C***********************************************************************
  45. -INC CCGEOME
  46. -INC PPARAM
  47. -INC CCOPTIO
  48. -INC CCREEL
  49. -INC SMCOORD
  50. -INC TMATOP1
  51. *-INC SMETRIQ
  52.       POINTEUR KCMETR.METRIQ
  53. *-INC SMELEMX
  54. -INC SMLREEL
  55.       PARAMETER(NMET=6)
  56.       DIMENSION XMET(NMET)
  57. *      DIMENSION XMET2(2,2)
  58. *      DIMENSION XMET3(3,3)
  59.       DIMENSION XMETV(3,3)
  60.       DIMENSION XJAC(3,3)
  61.       DIMENSION XJTMJ(3,3)
  62.       DIMENSION XMJ(3,3)
  63. *      REAL*8 DXI(3)
  64.       DIMENSION A(3,3),D(3)
  65. * Derivative of the affine barycentric map
  66.       DIMENSION DABM(4,3)
  67. * Simplex node coordinates
  68.       DIMENSION SNCO(3,4)
  69.  
  70.       DIMENSION IDXSYM(3,3,3)
  71. *
  72. * Statement functions
  73. *      DISTA(A,B,C,D)=SQRT((A-C)*(A-C)+(B-D)*(B-D))
  74. *      DISTB(A,B,C,D,E,F)=SQRT((A-D)*(A-D)+(B-E)*(B-E)+(C-F)*(C-F))
  75. *      DISTA(A,B)=SQRT(A*A+B*B)
  76. *      DISTB(A,B,C)=SQRT(A*A+B*B+C*C)
  77. *      DETTRI(A11,A12,A21,A22)=A11*A22-A21*A12
  78.       DETTET(A11,A12,A13,A21,A22,A23,A31,A32,A33)=
  79.      &  A11*(A22*A33-A23*A32)
  80.      &       +A12*(A23*A31-A21*A33)
  81.      &       +A13*(A21*A32-A22*A31)
  82. *
  83.       DATA ((XMETV(I,J),I=1,3),J=1,3) /9*0.D0/
  84.       DATA ((XJAC(I,J),I=1,3),J=1,3)  /9*0.D0/
  85.       DATA ((XJTMJ(I,J),I=1,3),J=1,3) /9*0.D0/
  86.       DATA ((XMJ(I,J),I=1,3),J=1,3) /9*0.D0/
  87.       DATA ((A(I,J),I=1,3),J=1,3) /9*0.D0/
  88.       DATA (D(I),I=1,3) /3*0.D0/
  89.       DATA ((DABM(I,J),I=1,4),J=1,3) /12*0.D0/
  90.       DATA ((SNCO(I,J),I=1,3),J=1,4) /12*0.D0/
  91.       DATA (((IDXSYM(I,J,K),I=1,1),J=1,1),K=1,1) /1/
  92.       DATA (((IDXSYM(I,J,K),I=1,2),J=1,2),K=2,2) /1,2,2,3/
  93.       DATA (((IDXSYM(I,J,K),I=1,3),J=1,3),K=3,3) /1,2,4,2,3,5,4,5,6/
  94. *
  95. *
  96. * Executable statements
  97. *
  98. *  Choix du critère sans metrique
  99. *  1 : 2D : D(2)/D(1), XALI
  100. *      3D : D(3)/D(1), D(2)/D(1), XALI
  101. *  2 : 2D : D(2)/D(1)
  102. *      3D : D(3)/D(1), D(2)/D(1)
  103. *  3 : 2D : XALI
  104. *      3D : XALI
  105. *  4 : 2D : sort(D(2)/D(1), XALI)
  106. *      3D : sort(D(3)/D(1), D(2)/D(1), XALI)
  107. *  5 : 2D : D(2)/D(1)
  108. *      3D : D(3)/D(1)
  109.       ICRIT=5
  110. * Choix du critere avec metrique
  111. * 1 : (SORT MIN(D(I),1./D(I))), XALI, D(IDIM)/D(1)
  112. * 2 : (SORT MIN(D(I),1./D(I))), XALI
  113. * 3 : (SORT MIN(D(I),1./D(I)), XALI)
  114. * 4 : (SORT MIN(D(I),1./D(I)), D(IDIM)/D(1))
  115. * 5 : (SORT MIN(D(I),1./D(I)), XALI, D(IDIM)/D(1))
  116. * 6 : (SORT MIN(D(I),1./D(I)))
  117. * 7 : (SORT MIN(D(I),1./D(I)), D(IDIM)/D(1)) si SORT < 0.5
  118. *     (SORT MIN(D(I),1./D(I)), XALI) si SORT > 0.5
  119. * 8 : (SORT_i-1 MIN(D(I),1./D(I)), D(IDIM)/D(1))
  120. * 9 : (SORT_i-1 MIN(D(I),1./D(I)), XALI)
  121. * 10 : (SORT_i-1 MIN(D(I),1./D(I))), D(IDIM)/D(1)
  122. * 11 : (SORT_i-1 MIN(D(I),1./D(I))), XALI
  123.       JCRIT=5
  124. *     write(ioimp,*) 'quali6: NKPVIR=',NKPVIR
  125.       IF (IMET.EQ.0) THEN
  126.          IF (ICRIT.EQ.1) THEN
  127.             ISTRID=IDIM
  128.          ELSEIF (ICRIT.EQ.2) THEN
  129.             ISTRID=MAX(1,IDIM-1)
  130.          ELSEIF (ICRIT.EQ.3.OR.ICRIT.EQ.5) THEN
  131.             ISTRID=1
  132.          ELSE
  133.             ISTRID=IDIM
  134.          ENDIF
  135.       ELSE
  136.          IF ((JCRIT.EQ.1).OR.(JCRIT.EQ.5)) THEN
  137.             ISTRID=IDIM+2
  138.          ELSEIF (JCRIT.EQ.6.OR.JCRIT.EQ.8.OR.JCRIT.EQ.9) THEN
  139.             ISTRID=IDIM
  140.          ELSE
  141.             ISTRID=IDIM+1
  142.          ENDIF
  143.       ENDIF
  144.       NDQC=0
  145.       IF (IELDEB.GT.IELFIN) RETURN
  146.       IDIMP1=IDIM+1
  147. *      NBNN=NUMX(/1)
  148.       NBNN=NNCOU
  149. *      NBELEM=NUMX(/2)
  150. *?      JG=PROG(/1)
  151. * Trop de verif nuit...
  152. *      IF (NBNN.NE.IDIMP1.OR.JG.NE.NBELEM) THEN
  153.       IF (NBNN.NE.IDIMP1) THEN
  154.          CALL ERREUR(5)
  155.          RETURN
  156.       ENDIF
  157.       IF
  158.      $     (.NOT.(IELDEB.GE.1.AND.IELFIN.GE.IELDEB.AND.NLCOU.GE.IELFIN
  159.      $     .AND.NUMX(/2).GE.NLCOU))  THEN
  160.          WRITE(IOIMP,*) 'coucou quali6'
  161.          write(ioimp,*) 'IELDEB=',IELDEB
  162.          write(ioimp,*) 'IELFIN=',IELFIN
  163.          write(ioimp,*) 'NLCOU=',NLCOU
  164.          write(ioimp,*) 'NUM2=',NUMX(/2)
  165.          write(ioimp,*) 'MELEMX=',MELEMX
  166.          CALL ECMELX(MELEMX,0)
  167.          CALL ERREUR(5)
  168.          RETURN
  169.       ENDIF
  170.  
  171. *     DO 10 IBELEM=1,NUMX(/2)
  172.       DO 10 IBELEM=IELDEB,IELFIN
  173. *     WRITE(IOIMP,*) 'IBELEM=',IBELEM
  174. *     Calcul de la métrique moyenne M soit dans XMETD (scalaire)
  175. *                                     soit dans XMETV (tenseur SPD)
  176. * Derivative of the affine barycentric map M : lambda(x) = Mx +c
  177. *     Les coordonnees barycentriques sont definies par rapport au
  178. *     simplex regulier de cote 1, centre sur l'origine. Le noeud sommet
  179. *     a toutes ses coordonnees nulles sauf la derniere
  180. * Initialisations au premier pas
  181.          IF (IBELEM.EQ.IELDEB) THEN
  182.             IF (IDIM.GE.1) THEN
  183.                DABM(1,1)=-1.D0
  184.                DABM(2,1)=+1.D0
  185.                IF (IDIM.GE.2) THEN
  186.                   DABM(1,2)=-1.D0/SQRT(3.D0)
  187.                   DABM(2,2)=-1.D0/SQRT(3.D0)
  188.                   DABM(3,2)=+2.D0/SQRT(3.D0)
  189.                   IF (IDIM.GE.3) THEN
  190.                      DABM(1,3)=-1.D0/SQRT(6.D0)
  191.                      DABM(2,3)=-1.D0/SQRT(6.D0)
  192.                      DABM(3,3)=-1.D0/SQRT(6.D0)
  193.                      DABM(4,3)=+3.D0/SQRT(6.D0)
  194.                   ENDIF
  195.                ENDIF
  196.             ENDIF
  197. *
  198.             IF (IMET.EQ.1) XMETD=1.D0/DENSIT
  199.             IF (IMET.EQ.2) XMETD=1.D0/XDENS
  200.             IF (IMET.EQ.3) NFMET=1
  201.             IF (IMET.EQ.4) NFMET=IDIM*(IDIM+1)/2
  202.          ENDIF
  203. *
  204.          if (imet.gt.0) then
  205.             IF (IMET.GE.1.AND.IMET.LE.3) THEN
  206.                IF (IMET.EQ.3) THEN
  207.                   YDENS=0.D0
  208.                   DO I=1,IDIMP1
  209.                      INO=NUMX(I,IBELEM)
  210.                      IF (NKPVIR.NE.0) THEN
  211.                         IF (INO.LE.NKPVIR) GOTO 10
  212.                      ENDIF
  213. *     Ici on fait la moyenne arithmétique
  214. *     mais kcmetr contient le log du tenseur si imomet=1
  215.                      YDENS=YDENS+KCMETR.XIN(1,INO)
  216.                   ENDDO
  217.                   YDENS=YDENS/IDIMP1
  218.                   if (imomet.eq.1) then
  219.                      YDENS=EXP(YDENS)
  220.                   endif
  221.                   XMETD2=YDENS
  222.                ELSE
  223.                   XMETD2=XMETD**2
  224.                ENDIF
  225.             ELSEIF (IMET.EQ.4) THEN
  226.                DO J=1,NFMET
  227.                   XMET(J)=0.D0
  228.                ENDDO
  229.                DO I=1,IDIMP1
  230.                   INO=NUMX(I,IBELEM)
  231.                   IF (NKPVIR.NE.0) THEN
  232.                      IF (INO.LE.NKPVIR) GOTO 10
  233.                   ENDIF
  234.                   DO J=1,NFMET
  235.                      XMET(J)=XMET(J)+KCMETR.XIN(J,INO)
  236.                   ENDDO
  237.                ENDDO
  238.                DO J=1,NFMET
  239.                   XMET(J)=XMET(J)/IDIMP1
  240.                ENDDO
  241. *
  242.                if (imomet.eq.1) then
  243.                   DO J=1,IDIM
  244.                      DO I=1,IDIM
  245.                         A(I,J)=XMET(IDXSYM(I,J,IDIM))
  246.                      ENDDO
  247.                   ENDDO
  248. *     Exponentielle du tenseur symétrique
  249.                   IOTENS=8
  250.                   IKAS=3
  251.                   CALL TENS2(IOTENS,IKAS,A,D,XMETV)
  252.                   IF (IERR.NE.0) RETURN
  253.                else
  254.                   DO J=1,IDIM
  255.                      DO I=1,IDIM
  256.                         XMETV(I,J)=XMET(IDXSYM(I,J,IDIM))
  257.                      ENDDO
  258.                   ENDDO
  259.                endif
  260.             ELSE
  261.                WRITE(IOIMP,*) 'quali6 imet=',IMET
  262.                CALL ERREUR(5)
  263.                RETURN
  264.             ENDIF
  265.          endif
  266. * Calcul du jacobien de la transformation geometrique entre
  267. * l'element regulier de coté 1 et l'element courant
  268. *   Coordonnees des noeuds
  269.          DO J=1,IDIMP1
  270.             INOD=NUMX(J,IBELEM)
  271.             IF (NKPVIR.NE.0) THEN
  272.                IF (INOD.LE.NKPVIR) goto 10
  273.             ENDIF
  274.             IPNOD=(INOD-1)*IDIMP1
  275.             DO I=1,IDIM
  276.                SNCO(I,J)=XCOOR(IPNOD+I)
  277.             ENDDO
  278.          ENDDO
  279. *         write(ioimp,*) 'SNCO,I,J=',IDIM,IDIMP1
  280. *         write(ioimp,*) ((SNCO(I,J),I=1,IDIM),J=1,IDIMP1)
  281. *         write(ioimp,*) 'DABM,I,J=',IDIMP1,IDIM
  282. *         write(ioimp,*) ((DABM(I,J),I=1,IDIMP1),J=1,IDIM)
  283. *   Matrice Jacobienne de la transformation J = SNCO*DABM
  284.          DO J=1,IDIM
  285.             DO I=1,IDIM
  286.                XIJ=0.D0
  287.                DO K=1,IDIMP1
  288.                   XIJ=XIJ+SNCO(I,K)*DABM(K,J)
  289.                ENDDO
  290.                XJAC(I,J)=XIJ
  291.             ENDDO
  292.          ENDDO
  293.          IF (IDIM.EQ.1) THEN
  294.             XDETJ=XJAC(1,1)
  295.          ELSEIF (IDIM.EQ.2) THEN
  296.             XDETJ=XJAC(1,1)*XJAC(2,2)-XJAC(2,1)*XJAC(1,2)
  297.          ELSEIF (IDIM.EQ.3) THEN
  298.             XDETJ=DETTET(XJAC(1,1),XJAC(1,2),XJAC(1,3),XJAC(2
  299.      $           ,1),XJAC(2,2),XJAC(2,3),XJAC(3,1),XJAC(3,2)
  300.      $           ,XJAC(3,3))
  301.          ELSE
  302.             INTERR(1)=IDIM
  303.             CALL ERREUR(709)
  304.             RETURN
  305.          ENDIF
  306. *         write(ioimp,*) 'XJAC,I,J=',IDIM,IDIM
  307. *         write(ioimp,*) ((XJAC(I,J),I=1,IDIM),J=1,IDIM)
  308. * Matrice JtMJ
  309.          IF (IMET.LT.4) THEN
  310.             DO K=1,IDIM
  311.                DO I=1,IDIM
  312.                   XIK=0.D0
  313.                   DO J=1,IDIM
  314.                      XIK=XIK+XJAC(J,I)*XJAC(J,K)
  315.                   ENDDO
  316.                   IF (IMET.EQ.0) THEN
  317.                      XJTMJ(I,K)=XIK
  318.                   ELSE
  319.                      XJTMJ(I,K)=XIK*XMETD2
  320.                   ENDIF
  321.                ENDDO
  322.             ENDDO
  323.             IF (IMET.EQ.0) THEN
  324.                XDETM=1.D0
  325.             ELSE
  326.                XDETM=XMETD2**IDIM
  327.             ENDIF
  328.          ELSE
  329.             DO L=1,IDIM
  330.                DO J=1,IDIM
  331.                   XJL=0.D0
  332.                   DO K=1,IDIM
  333. * Utilisons la symetrie de XMETV
  334.                      XJL=XJL+XMETV(K,J)*XJAC(K,L)
  335.                   ENDDO
  336.                   XMJ(J,L)=XJL
  337.                ENDDO
  338.             ENDDO
  339.             DO L=1,IDIM
  340.                DO I=1,IDIM
  341.                   XIL=0.D0
  342.                   DO J=1,IDIM
  343.                      XIL=XIL+XJAC(J,I)*XMJ(J,L)
  344.                   ENDDO
  345.                   XJTMJ(I,L)=XIL
  346.                ENDDO
  347.             ENDDO
  348.             IF (IDIM.EQ.1) THEN
  349.                XDETM=XMETV(1,1)
  350.             ELSEIF (IDIM.EQ.2) THEN
  351.                XDETM=XMETV(1,1)*XMETV(2,2)-XMETV(2,1)*XMETV(1,2)
  352.             ELSEIF (IDIM.EQ.3) THEN
  353.                XDETM=DETTET(XMETV(1,1),XMETV(1,2),XMETV(1,3),XMETV(2
  354.      $              ,1),XMETV(2,2),XMETV(2,3),XMETV(3,1),XMETV(3,2)
  355.      $              ,XMETV(3,3))
  356.             ELSE
  357.                INTERR(1)=IDIM
  358.                CALL ERREUR(709)
  359.                RETURN
  360.             ENDIF
  361.          ENDIF
  362. *         WRITE(IOIMP,*) 'XDETM=',XDETM
  363. *         write(ioimp,*) 'XJTMJ,I,J=',IDIM,IDIM
  364. *     write(ioimp,*) ((XJTMJ(I,J),I=1,IDIM),J=1,IDIM)
  365.          XDETJM=XDETJ*SQRT(XDETM)
  366. * Determinant et trace de JTMJ
  367.          IF (IDIM.EQ.1) THEN
  368.             D(1)=XDETJM**2
  369.          ELSEIF (IDIM.EQ.2) THEN
  370. *            XDET=XJTMJ(1,1)*XJTMJ(2,2)-XJTMJ(2,1)*XJTMJ(1,2)
  371.             CALL JACOD2(XJTMJ,D)
  372.          ELSEIF (IDIM.EQ.3) THEN
  373. *            XDET=DETTET(XJTMJ(1,1),XJTMJ(1,2),XJTMJ(1,3),XJTMJ(2
  374. *     $           ,1),XJTMJ(2,2),XJTMJ(2,3),XJTMJ(3,1),XJTMJ(3,2)
  375. *     $           ,XJTMJ(3,3))
  376.             CALL JACOD3(XJTMJ,3,D)
  377.          ELSE
  378.             WRITE(IOIMP,*) 'quali6 idim=',IDIM
  379.             INTERR(1)=IDIM
  380.             CALL ERREUR(709)
  381.             RETURN
  382.          ENDIF
  383.          XPET=XPETIT*10.D0
  384.          XTR=D(1)
  385.          IF (IDIM.EQ.1) THEN
  386.             XALIN1=1
  387.          ELSE
  388.             DO I=2,IDIM
  389.                XTR=XTR+D(I)
  390.             ENDDO
  391.             XL2T=XTR/IDIM
  392.             IF (IDIM.EQ.2) THEN
  393.                XLTD=XL2T
  394.             ELSEIF (IDIM.EQ.3) THEN
  395.                XLTD=XL2T*SQRT(XL2T)
  396.             ELSE
  397.                INTERR(1)=IDIM
  398.                CALL ERREUR(709)
  399.                RETURN
  400.             ENDIF
  401.             XALIN1=ABS(XDETJM)/(MAX(XLTD,XPET))
  402. *     XALIN1=sqrt(max(D(IDIM)/(MAX(D(1),XPET)),xzero))
  403.             IF (IDIM.EQ.3) XALIN1=SQRT(XALIN1)
  404.          ENDIF
  405. * Par rapport au livre de Huang p.205, XALIN1 vaut 1 / (Qali^(n-1)) (n>=2)
  406. * Comme Qali est minore par le rapport d'aspect d'un element, il
  407. * faut comparer XALIN1 a des (rapports de longueur)^(n-1)
  408. * Il y a donc un carre en dimension 3 cf. les expressions de XQUALN
  409. * plus bas (voir aussi deadutil.procedur qui exprime des indicateurs
  410. * en rapports de longueur)
  411. * On a choisi d'exprimer XALIN1 en fonction de XDETJM directement
  412. * car la presque nullite de XDETJ permet de detecter les elements
  413. * plats.
  414. * Si on l'eleve a la puissance (1/3), ca ne marche pas.
  415. * On devra sans doute faire mieux pour etre vraiment robuste....
  416. *         write(ioimp,*) 'XALIN1=',XALIN1
  417. * Les valeurs propres sont censees etre positives mais pas garanti
  418. * donc on prend la valeur absolue et on reclasse
  419.          JELDEB=(IELDEB-1+NDQC)*ISTRID
  420. *         WRITE(IOIMP,*) '1',(D(II),II=1,IDIM)
  421.          DO I=1,IDIM
  422. *           D(I)=ABS(D(I))
  423.             D(I)=SQRT(MAX(D(I),XZERO))
  424.          ENDDO
  425.          IF (IMET.EQ.0) THEN
  426.             IF (IDIM.EQ.1) THEN
  427.                PROG(JELDEB+1)=1.D0
  428.             ELSE
  429.                IF (ICRIT.NE.3) THEN
  430.                   IF (ICRIT.EQ.5) THEN
  431.                      KMAX=1
  432.                   ELSE
  433.                      KMAX=IDIM-1
  434.                   ENDIF
  435.                   IF (D(1).NE.XZERO) THEN
  436. *                D(1)=MAX(D(1),XPET)
  437.                      DO K=1,KMAX
  438.                         PROG(JELDEB+K)=D(IDIM-K+1)/D(1)
  439.                      ENDDO
  440.                   ELSE
  441.                      DO K=1,KMAX
  442.                         PROG(JELDEB+K)=XZERO
  443.                      ENDDO
  444.                   ENDIF
  445.                ENDIF
  446.                IF (ICRIT.NE.2) THEN
  447.                   PROG(JELDEB+ISTRID)=XALIN1
  448.                ENDIF
  449.                IF (ICRIT.EQ.4) THEN
  450.                   CALL ORDO01(PROG(JELDEB+1),IDIM,.TRUE.)
  451.                ENDIF
  452.             ENDIF
  453.          ELSE
  454.             IF (JCRIT.EQ.1.OR.JCRIT.EQ.4.OR.JCRIT.EQ.5) THEN
  455.                IF (D(1).NE.XZERO) THEN
  456. *               PROG(JELDEB+ISTRID)=D(IDIM)/MAX(D(1),XPET)
  457.                   PROG(JELDEB+ISTRID)=D(IDIM)/D(1)
  458.                ELSE
  459.                   PROG(JELDEB+ISTRID)=XZERO
  460.                ENDIF
  461.             ENDIF
  462. *            WRITE(IOIMP,*) '3',(D(II),II=1,IDIM)
  463.             DO K=1,IDIM
  464.                XL=D(K)
  465.                XLI=1.D0/MAX(D(K),XPET)
  466.                D(K)=MIN(XL,XLI)
  467. * Plante avec le compilo du jour
  468. *     IF (XL.GT.1.D0) THEN
  469. **                  WRITE(IOIMP,*) 'XL=',XL
  470. *                  D(K)=1.D0/XL
  471. *               ENDIF
  472.             ENDDO
  473.             DO K=1,IDIM
  474.                PROG(JELDEB+K)=D(K)
  475.             ENDDO
  476.             IF ((JCRIT.EQ.1).OR.(JCRIT.EQ.5)) THEN
  477.                PROG(JELDEB+ISTRID-1)=XALIN1
  478.             ELSEIF ((JCRIT.EQ.2).OR.(JCRIT.EQ.3)) THEN
  479.                PROG(JELDEB+ISTRID)=XALIN1
  480.             ENDIF
  481.             IF ((JCRIT.EQ.1).OR.(JCRIT.EQ.2).OR.(JCRIT.GE.6)) THEN
  482.                CALL ORDO01(PROG(JELDEB+1),IDIM,.TRUE.)
  483.             ELSEIF ((JCRIT.EQ.3).OR.(JCRIT.EQ.4)) THEN
  484.                CALL ORDO01(PROG(JELDEB+1),IDIM+1,.TRUE.)
  485.             ELSEIF (JCRIT.EQ.5) THEN
  486.                CALL ORDO01(PROG(JELDEB+1),IDIM+2,.TRUE.)
  487.             ENDIF
  488.             IF (JCRIT.EQ.7) THEN
  489.                IF (PROG(JELDEB+1).LT.1.05D0) THEN
  490.                   IF (D(1).NE.XZERO) THEN
  491.                      PROG(JELDEB+ISTRID)=D(IDIM)/D(1)
  492.                   ELSE
  493.                      PROG(JELDEB+ISTRID)=XZERO
  494.                   ENDIF
  495.                ELSE
  496.                   PROG(JELDEB+ISTRID)=XALIN1
  497.                ENDIF
  498.                CALL ORDO01(PROG(JELDEB+1),IDIM+1,.TRUE.)
  499.             ENDIF
  500.             IF (JCRIT.EQ.8.OR.JCRIT.EQ.10) THEN
  501.                IF (D(1).NE.XZERO) THEN
  502.                   PROG(JELDEB+ISTRID)=D(IDIM)/D(1)
  503.                ELSE
  504.                   PROG(JELDEB+ISTRID)=XZERO
  505.                ENDIF
  506.                IF (JCRIT.EQ.8)
  507.      $              CALL ORDO01(PROG(JELDEB+1),IDIM,.TRUE.)
  508.             ENDIF
  509.             IF (JCRIT.EQ.9.OR.JCRIT.EQ.11) THEN
  510.                PROG(JELDEB+ISTRID)=XALIN1
  511.                IF (JCRIT.EQ.9)
  512.      $              CALL ORDO01(PROG(JELDEB+1),IDIM,.TRUE.)
  513.             ENDIF
  514.          ENDIF
  515. * Scaling ? mmmmm, petit doute
  516.          NDQC=NDQC+1
  517.  10   CONTINUE
  518.       RETURN
  519. *
  520. * formats
  521. *
  522.   188  FORMAT (2X,12(A6,'=',1PG12.5,2X))
  523. *
  524. * End of subroutine QUALI6
  525. *
  526.       END
  527.  
  528.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales