Télécharger la plaquette Cast3M


fr en

Télécharger yjohns.eso

Retour à la liste

Numérotation des lignes : 

yjohns

  1. C YJOHNS    SOURCE    CB215821  26/08/24    21:19:00     12622          
  2.       SUBROUTINE YJOHNS
  3. C
  4. C     VERSION VECTORISEE
  5. C
  6. C     Les éléments sont groupés en paquets de LRV éléments, LRV étant
  7. C     la longueur des registres vectoriels de la machine cible, i.e
  8. C     64 sur Cray, 128 ou 256 sur IBM 3090VF. On promène une fenêtre
  9. C     de longueur LRV sur la boucle générale de longueur NEL.
  10. C
  11.  
  12. C    KDESIGN n'est pas défini  !!!!!!
  13. C
  14. C
  15.      & (HR,RPG,DRR,LE,NEL,K0,NPTD,IES,NP,IAXI,
  16.      & IPADI,IKOMP,IKAS,
  17.      & ALFE,IND1,UN,INDU,NPTS,IPADS,
  18.      & TN,QE,IKS,
  19.      & HRN,G,ALT,SGT,
  20.      & VOLU,COTE,NELZ,
  21.      & DTM1,DT,DTT1,DTT2,NUEL,DIAEL)
  22.  
  23.       IMPLICIT INTEGER(I-N)
  24.       IMPLICIT REAL*8 (A-H,O-Z)
  25.  
  26. C***********************************************************************
  27. C
  28. C  CE SP DISCRETISE UNE EQUATION GENERALE DE TRANSPORT-DIFFUSION AVEC
  29. C  SOURCE.
  30. C    EN 2D SUR LES ELEMENTS QUA4 ET TRI3      PLAN OU AXI
  31. C    EN 3D SUR LES ELEMENTS CUB8 ET PRI6
  32. C  LES OPERATEURS SONT "SOUS-INTEGRES"
  33. C
  34. C
  35. C                        APPELE PAR YTSCAL
  36. C
  37. C
  38. C***********************************************************************
  39.  
  40. -INC CCVQUA4
  41. -INC CCREEL
  42. -INC SMCOORD
  43.  
  44. -INC PPARAM
  45. -INC CCOPTIO
  46.  
  47. C
  48. C Longueur des registres vectoriels de la machine cible
  49. C On prend 64 pour ne pas augmenter la taille des tableaux
  50. C nécessaires à la vectorisation.
  51. C
  52.       PARAMETER(LRV=64)
  53.  
  54.       DIMENSION UN(NPTS,IES),HRN(NPTD),TN(NPTD)
  55.       DIMENSION COTE(NELZ,IES),VOLU(NELZ),QE(*)
  56.       DIMENSION ALFE(*),ALT(*)
  57.  
  58.       DIMENSION IPADI(1),LE(NP,1),IPADS(*)
  59.       DIMENSION HR(NEL,NP,IES),RPG(1),DRR(NP,NEL)
  60.  
  61.       DIMENSION BF(9,9)
  62.       DIMENSION QGGT(8,8),Q1(8,8),Q2(8,8),Q3(8,8)
  63.  
  64.       DIMENSION COEF(LRV),AIRE(LRV)
  65.       DIMENSION AL2(LRV),AH2(LRV),AP2(LRV)
  66.       DIMENSION AL(LRV),AH(LRV),AP(LRV)
  67.       DIMENSION ALFT(LRV),QT(LRV)
  68.       DIMENSION XMB(LRV),XMH(LRV)
  69.       DIMENSION CFM(LRV)
  70. C     DIMENSION CF1(LRV),CF2(LRV),CF3(LRV)
  71.       DIMENSION DR(LRV,9)
  72.       DIMENSION UM(LRV),UP(LRV)
  73.       DIMENSION UIX(LRV,9),UIY(LRV,9),UIZ(LRV,9)
  74.       DIMENSION TETAC(LRV,9),TETAD(LRV,9),TETA(LRV,9)
  75.       DIMENSION UMI(LRV,3),UPI(LRV,3)
  76.       DIMENSION CXT(LRV),CYT(LRV),CXY(LRV)
  77.       DIMENSION DXT(LRV),DYT(LRV),DXY(LRV)
  78.       DIMENSION CZT(LRV),CXZ(LRV),CYZ(LRV)
  79.       DIMENSION DZT(LRV),DXZ(LRV),DYZ(LRV)
  80.       DIMENSION SBF(LRV,9)
  81.       DIMENSION GRADT(LRV,3)
  82.       DIMENSION BM(LRV),BP(LRV)
  83.  
  84.       DIMENSION BMX(LRV),BMY(LRV),BPX(LRV),BPY(LRV)
  85.  
  86.       REAL*8    G(*),n(2,4),MMAX(4),lll,Fi
  87.       REAL*8    b(2),Nm,x(4),y(4),kkk(4),Kchap
  88.       INTEGER   zz,p,kdesign
  89.       PARAMETER (zz=1)
  90.  
  91.       SAVE IPAS,QGGT,Q1,Q2,Q3
  92.  
  93.       DATA CD/1.D0/
  94.  
  95.       DATA IPAS/0/
  96.  
  97.       DATA IDCENN/2/
  98. C************************************************************************
  99. C
  100. C           INITIALISATIONS DIVERSES
  101. C
  102.       KDESIGN=0
  103.       NK=K0
  104.  
  105. C                           ********
  106. C                           *  2D  *
  107. C                           ********
  108.  
  109.  
  110.       IF (NP.EQ.3) THEN
  111.       x(4)=0.d0
  112.       y(4)=0.d0
  113.       ENDIF
  114.  
  115.       IF(IES.EQ.3)GO TO 10
  116.  
  117.       IAX1=0
  118.       IAX2=0
  119.       IF(IAXI.EQ.1)IAX2=1
  120.       IF(IAXI.EQ.2)IAX1=1
  121.  
  122.       HMIN = 1.D20
  123.       HMAX = 0.D0
  124.  
  125. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  126.  
  127. C                           DIFFERENCES TRIANGLE / QUADRANGLE
  128.       QUA4=0.D0
  129.       IF(NP.EQ.4)QUA4=1.D0
  130. C
  131. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  132.  
  133. C
  134. C Calcul du nombre de paquets de LRV éléments
  135. C
  136.       NNN=MOD(NEL,LRV)
  137.       IF(NNN.EQ.0) NPACK=NEL/LRV
  138.       IF(NNN.NE.0) NPACK=1+(NEL-NNN)/LRV
  139.       KPACKD=1
  140.       KPACKF=NPACK
  141. C
  142. C ******* BOUCLE SUR LES PAQUETS DE LRV ELEMENTS **********
  143. C
  144.       DO 7001 KPACK=KPACKD,KPACKF
  145. C
  146. C ======= A L'INTERIEUR DE CHAQUE PAQUET DE LRV ELEMENTS =======
  147. C
  148. C 1. Calcul des limites du paquet courant.
  149.       KDEB=1+(KPACK-1)*LRV
  150.       KFIN=MIN(NEL,KDEB+LRV-1)
  151. C
  152.       DO 7002 K=KDEB,KFIN
  153.       KP=K-KDEB+1
  154.       NK=K+K0
  155.       NK1=(1-IND1)*(NK-1)+1
  156.       ALFT(KP)=ALFE(NK1)+XPETIT
  157.       AIRE(KP)=VOLU(NK)
  158.       AL(KP)=COTE(NK,1)
  159.       AH(KP)=COTE(NK,2)
  160.       AL2(KP)=1.D0/AL(KP)/AL(KP)
  161.       AH2(KP)=1.D0/AH(KP)/AH(KP)
  162.       XMH(KP)=(AL(KP)+AH(KP))/2.D0
  163.  7002 CONTINUE
  164.  
  165.       IF((IKOMP.EQ.0.AND.IKAS.EQ.5).OR.
  166.      &(IKOMP.EQ.1.AND.IKAS.EQ.6))THEN
  167.       DO 7003 K=KDEB,KFIN
  168.       KP=K-KDEB+1
  169.       NK=K+K0
  170.       ALFT(KP)=ALFT(KP)+ALT(NK)/SGT
  171.  7003 CONTINUE
  172.       ENDIF
  173.  
  174. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  175.  
  176. C Initialisation des UMI avant accumulation
  177.       DO 7005 K=KDEB,KFIN
  178.       KP=K-KDEB+1
  179.       UMI(KP,1)=XPETIT
  180.       UMI(KP,2)=XPETIT
  181.       UPI(KP,1)=XPETIT
  182.       UPI(KP,2)=XPETIT
  183.       GRADT(KP,1)=XPETIT
  184.       GRADT(KP,2)=XPETIT
  185.  7005 CONTINUE
  186.  
  187.       DO 81082 I=1,NP
  188. C*IBMDIR* PREFER VECTOR
  189.       DO 7006 K=KDEB,KFIN
  190.       KP=K-KDEB+1
  191.       NF=IPADI(LE(I,K))
  192.       NFU=IPADS(LE(I,K))
  193.       NFU=(1-INDU)*(NFU-1)+1
  194.       DR(KP,I)=DRR(I,K)
  195.       UIX(KP,I)=UN(NFU,1)
  196.       UIY(KP,I)=UN(NFU,2)
  197.       UMI(KP,1)=UMI(KP,1)+UN(NFU,1)*DR(KP,I)
  198.       UMI(KP,2)=UMI(KP,2)+UN(NFU,2)*DR(KP,I)
  199.  
  200.       TETAC(KP,I)=HRN(NF)
  201.       TETAD(KP,I)=TN(NF)
  202.       GRADT(KP,1)=GRADT(KP,1)+HR(K,I,1)*TETAC(KP,I)
  203.       GRADT(KP,2)=GRADT(KP,2)+HR(K,I,2)*TETAC(KP,I)
  204.  7006 CONTINUE
  205. 81082 CONTINUE
  206.  
  207. C
  208. C Initialisation de la variable d'accumulation SBF au terme source
  209. C
  210. C       write(6,*)' IKomp,ias=',IKomp,ikas
  211. C       write(6,*)' IKS,IND1,INDU=',IKS,IND1,INDU
  212.  
  213.       IF(IKOMP.EQ.0)THEN
  214.  
  215.          DO 70021 K=KDEB,KFIN
  216.          KP=K-KDEB+1
  217.          NK=K+K0
  218.          NKS=(1-IKS)*(NK-1)+1
  219.          QT(KP)=QE(NKS)
  220. 70021    CONTINUE
  221.  
  222.          DO 81083 I=1,NP
  223. C*IBMDIR* PREFER VECTOR
  224.          DO 70062 K=KDEB,KFIN
  225.          KP=K-KDEB+1
  226.          SBF(KP,I)=-QT(KP)*DR(KP,I)
  227. 70062    CONTINUE
  228. 81083    CONTINUE
  229.  
  230.       ELSEIF(IKOMP.EQ.1)THEN
  231.  
  232.          DO 70023 K=KDEB,KFIN
  233.          KP=K-KDEB+1
  234.          NK=K+K0
  235.          NKS=(1-IKS)*(NK-1)+1
  236.          QT(KP)=QE(NKS)
  237. 70023    CONTINUE
  238.  
  239.          DO 81084 I=1,NP
  240. C*IBMDIR* PREFER VECTOR
  241.          DO 70064 K=KDEB,KFIN
  242.          KP=K-KDEB+1
  243.          SBF(KP,I)=-QT(KP)*DR(KP,I)
  244. 70064    CONTINUE
  245. 81084    CONTINUE
  246.  
  247.       ENDIF
  248.  
  249.       DO 7007 K=KDEB,KFIN
  250.       KP=K-KDEB+1
  251.       UMI(KP,1)=UMI(KP,1)/AIRE(KP)
  252.       UMI(KP,2)=UMI(KP,2)/AIRE(KP)
  253.       UM(KP)=UMI(KP,1)*UMI(KP,1)+UMI(KP,2)*UMI(KP,2)
  254.       UM(KP)=SQRT(UM(KP)) + XPETIT
  255.  
  256.       A=GRADT(KP,1)*GRADT(KP,1)+GRADT(KP,2)*GRADT(KP,2)
  257.       A=SQRT(A)+XPETIT
  258.       GRADT(KP,1)=GRADT(KP,1)/A
  259.       GRADT(KP,2)=GRADT(KP,2)/A
  260.       UP(KP)=GRADT(KP,1)*UMI(KP,1)+GRADT(KP,2)*UMI(KP,2)
  261.       UPI(KP,1)=UP(KP)*GRADT(KP,1)
  262.       UPI(KP,2)=UP(KP)*GRADT(KP,2)
  263.       UP(KP) = ABS(UP(KP)) + XPETIT
  264.  7007 CONTINUE
  265.  
  266. *      DO 70074 K=KDEB,KFIN
  267. *      KP=K-KDEB+1
  268. *      BMX=UMI(KP,1)/AL(KP)
  269. *      BMY=UMI(KP,2)/AH(KP)
  270. *      BM(KP)=BMX*BMX+BMY*BMY
  271. *      BM(KP)=SQRT(BM(KP))+XPETIT
  272. *      BPX=UPI(KP,1)/AL(KP)
  273. *      BPY=UPI(KP,2)/AH(KP)
  274. *      BP(KP)=BPX*BPX+BPY*BPY
  275. *      BP(KP)=SQRT(BP(KP))+XPETIT
  276. *70074 CONTINUE
  277.  
  278.  
  279.       DO 70074 K=KDEB,KFIN
  280.       KP=K-KDEB+1
  281.       BMX(KP)=UMI(KP,1)/AL(KP)
  282.       BMY(KP)=UMI(KP,2)/AH(KP)
  283.  
  284.       BM(KP)=BMX(KP)*BMX(KP)+BMY(KP)*BMY(KP)
  285.       BM(KP)=SQRT(BM(KP))+XPETIT
  286.       BPX(KP)=UPI(KP,1)/AL(KP)
  287.       BPY(KP)=UPI(KP,2)/AH(KP)
  288.       BP(KP)=BPX(KP)*BPX(KP)+BPY(KP)*BPY(KP)
  289.       BP(KP)=SQRT(BP(KP))+XPETIT
  290. 70074 CONTINUE
  291.  
  292.  
  293. C
  294. C     AVANT DECENTREMENT
  295. C
  296.  
  297.  
  298.       IF(IKOMP.EQ.0)THEN
  299.  
  300. C ALFA,AKSI,CCU,CCT utilisées seulement dans la boucle 7008
  301.       IF(IDCENN.EQ.1)THEN
  302.       DO 70081 K=KDEB,KFIN
  303.       KP=K-KDEB+1
  304.       XMB(KP)=XMH(KP)
  305.       CXT(KP)=0.D0
  306.       CYT(KP)=0.D0
  307.       CXY(KP)=0.D0
  308. 70081 CONTINUE
  309.  
  310.       ELSEIF(IDCENN.EQ.2)THEN
  311.  
  312.       DO 7008 K=KDEB,KFIN
  313.       KP=K-KDEB+1
  314.  
  315.       HMK=2.D0*UM(KP)/BM(KP)
  316.       HPK=2.D0*UP(KP)/BP(KP)
  317.  
  318.       IF (HMK.LT.HMIN) HMIN=HMK
  319.       IF (HMK.GT.HMAX) HMAX=HMK
  320.  
  321.       XMB(KP)=HMK
  322.       ALFA=UM(KP)*HMK/ALFT(KP)/2.D0
  323.       AKSI=ALFA/3.D0
  324.       IF(ALFA.GT.3.D0)AKSI=1.D0
  325.       CCT=AKSI*HMK/2.d0/UM(KP)
  326.  
  327.       Fi= (GRADT(KP,1)*UMI(KP,1)
  328.      &+    GRADT(KP,2)*UMI(KP,2))
  329.  
  330.       Kchap=0.5*(XMB(KP)*abs(Fi))/(sqrt((GRADT(KP,1)**2+
  331.      &     GRADT(KP,2)**2))+XMB(KP))
  332.  
  333.       CXT(KP)=(UMI(KP,1)*UMI(KP,1)*CCT+Kchap)
  334.       CYT(KP)=(UMI(KP,2)*UMI(KP,2)*CCT+Kchap)
  335.       CXY(KP)=(UMI(KP,1)*UMI(KP,2)*CCT)
  336.  
  337. 7008  CONTINUE
  338.  
  339.       ELSEIF(IDCENN.EQ.3)THEN
  340.       DO 7018 K=KDEB,KFIN
  341.       KP=K-KDEB+1
  342. C DESIGN TEZDUYAR
  343.  
  344.       IF (KDESIGN.EQ.1) THEN
  345. *      IF (NP.NE.3) RETURN
  346.       S=0.d0
  347.       do 119 I=1,NP
  348.       S=S+abs(UMI(KP,1)*HR(K,I,1)+UMI(KP,2)*HR(K,I,2))/UM(KP)
  349.   119 CONTINUE
  350.       HMK=2.d0/S
  351. C DESIGN ORIGINAL
  352.       ELSEIF (KDESIGN.EQ.0) THEN
  353.       HMK=2.D0*UM(KP)/BM(KP)
  354.  
  355.  
  356.       ELSEIF (KDESIGN.EQ.2) THEN
  357.       do 774 i=1,NP
  358.       x(i)=XCOOR((LE(i,K)-1)*(IES+1)+1)
  359.       y(i)=XCOOR((LE(i,K)-1)*(IES+1)+2)
  360.   774 CONTINUE
  361.  
  362.       DJ=(x(2)-x(1)+(x(3)-x(4))*(NP-3))*(y(3)-y(1)+(y(4)-y(2))*(NP-3))
  363.      *  -(y(2)-y(1)+(y(3)-y(4))*(NP-3))*(x(3)-x(1)+(x(4)-x(2))*(NP-3))
  364.  
  365.       IF (NP.EQ.3) lll=1.d0
  366.       IF (NP.EQ.4) lll=4.d0
  367.       DJ=DJ/(lll**2)
  368.  
  369.       b(1)=(UMI(KP,1)*(y(3)-y(1)+(y(4)-y(2))*(NP-3))-UMI(KP,2)*(x(3)-
  370.      * x(1)+(x(4)-x(2))*(NP-3)))/DJ/lll
  371.       b(2)=(-UMI(KP,1)*(y(2)-y(1)+(y(3)-y(4))*(NP-3))+UMI(KP,2)*(x(2)-
  372.      * x(1)+(x(3)-x(4))*(NP-3)))/DJ/lll
  373.  
  374.       p=1
  375.       s=0.d0
  376.       do 304 kk=1,2
  377.       s=s+(abs(b(kk)))**p
  378.   304 CONTINUE
  379.       Nm=s**(1.d0/p)
  380.       HMK=2.d0*UM(KP)/Nm
  381.  
  382.       ELSEIF (KDESIGN.EQ.3) THEN
  383.       do 775 i=1,NP
  384.       x(i)=XCOOR((LE(i,K)-1)*(IES+1)+1)
  385.       y(i)=XCOOR((LE(i,K)-1)*(IES+1)+2)
  386.   775 CONTINUE
  387.  
  388.       DJ=(x(2)-x(1)+(x(3)-x(4))*(NP-3))*(y(3)-y(1)+(y(4)-y(2))*(NP-3))
  389.      *  -(y(2)-y(1)+(y(3)-y(4))*(NP-3))*(x(3)-x(1)+(x(4)-x(2))*(NP-3))
  390.  
  391.       IF (NP.EQ.3) lll=1.d0
  392.       IF (NP.EQ.4) lll=4.d0
  393.       DJ=DJ/(lll**2)
  394.  
  395.       b(1)=(UMI(KP,1)*(y(3)-y(1)+(y(4)-y(2))*(NP-3))-UMI(KP,2)*(x(3)-
  396.      * x(1)+(x(4)-x(2))*(NP-3)))/DJ/lll
  397.       b(2)=(-UMI(KP,1)*(y(2)-y(1)+(y(3)-y(4))*(NP-3))+UMI(KP,2)*(x(2)-
  398.      * x(1)+(x(3)-x(4))*(NP-3)))/DJ/lll
  399.  
  400.       p=2
  401.       s=0.d0
  402.       do 305 kk=1,2
  403.       s=s+(abs(b(kk)))**p
  404.   305 CONTINUE
  405.       Nm=s**(1.d0/p)
  406.       HMK=2.d0*UM(KP)/Nm
  407.  
  408.       ELSEIF (KDESIGN.EQ.4) THEN
  409.       do 776 i=1,NP
  410.       x(i)=XCOOR((LE(i,K)-1)*(IES+1)+1)
  411.       y(i)=XCOOR((LE(i,K)-1)*(IES+1)+2)
  412.   776 CONTINUE
  413.  
  414.       do 318 i=1,NP
  415.       kkk(i)=Aire(KP)*(UMI(KP,1)*HR(K,i,1)+UMI(KP,2)*HR(K,i,2))
  416.   318 CONTINUE
  417.  
  418.       s0=0.d0
  419.       s1=0.d0
  420.       do 319 i=1,NP
  421.       s0=s0+min(0.d0,kkk(i))
  422.       s1=s1+max(0.d0,kkk(i))
  423.   319 CONTINUE
  424.  
  425.       s2=0.d0
  426.       s3=0.d0
  427.       do 320 i=1,NP
  428.       s2=s2-min(0.d0,kkk(i))*x(i)/s0+max(0.d0,kkk(i))*x(i)/s1
  429.       s3=s3-min(0.d0,kkk(i))*y(i)/s0+max(0.d0,kkk(i))*y(i)/s1
  430.   320 CONTINUE
  431.       HMK=sqrt(s2**2+s3**2)
  432.  
  433.       ELSEIF (KDESIGN.EQ.5) THEN
  434.       do 800 i=1,NP
  435.       n(1,i)=abs(2.d0*Aire(KP)*HR(K,i,1))
  436.       n(2,i)=abs(2.d0*Aire(KP)*HR(K,i,2))
  437.       MMAX(i)=n(1,i)
  438.       IF (n(1,i).lt.n(2,i)) MMAX(i)=n(2,i)
  439.   800 CONTINUE
  440.  
  441.       HMK=0.d0
  442.       do 801 i=1,NP
  443.       IF (MMAX(i).gt.HMK) HMK=MMAX(i)
  444.   801 CONTINUE
  445.  
  446.       ENDIF
  447.  
  448.  
  449.       IF (HMK.LT.HMIN) HMIN=HMK
  450.       IF (HMK.GT.HMAX) HMAX=HMK
  451.  
  452.       XMB(KP)=HMK
  453.       ALFA=UM(KP)*HMK/ALFT(KP)/2.D0
  454.       AKSI=ALFA/3.D0
  455.       IF(ALFA.GT.3.D0)AKSI=1.D0
  456.       CCT=AKSI*HMK/2.d0/UM(KP)
  457.  
  458.       CXT(KP)=UMI(KP,1)*UMI(KP,1)*CCT
  459.       CYT(KP)=UMI(KP,2)*UMI(KP,2)*CCT
  460.       CXY(KP)=UMI(KP,1)*UMI(KP,2)*CCT
  461.  
  462. 7018  CONTINUE
  463.  
  464.       ELSEIF(IDCENN.EQ.4)THEN
  465.       DT19=DTM1*0.5D0
  466.       DO 7009 K=KDEB,KFIN
  467.       KP=K-KDEB+1
  468.       XMB(KP)=XMH(KP)
  469.       CXT(KP)=UMI(KP,1)*UMI(KP,1)*DT19
  470.       CYT(KP)=UMI(KP,2)*UMI(KP,2)*DT19
  471.       CXY(KP)=UMI(KP,1)*UMI(KP,2)*DT19
  472. 7009  CONTINUE
  473.  
  474.  
  475.       ENDIF
  476.       ELSEIF(IKOMP.EQ.1)THEN
  477.  
  478. C ALFA,AKSI,CCU,CCT utilisées seulement dans la boucle 7008
  479.       IF(IDCENN.EQ.1)THEN
  480.       DO 71081 K=KDEB,KFIN
  481.       KP=K-KDEB+1
  482.       XMB(KP)=XMH(KP)
  483.  
  484.       CXT(KP)=0.D0
  485.       CYT(KP)=0.D0
  486.       CXY(KP)=0.D0
  487.  
  488.       DXT(KP)=0.D0
  489.       DYT(KP)=0.D0
  490.       DXY(KP)=0.D0
  491.  
  492. 71081 CONTINUE
  493.  
  494.       ELSEIF(IDCENN.EQ.2)THEN
  495.       DO 7108 K=KDEB,KFIN
  496.       KP=K-KDEB+1
  497. C DESIGN TEZDUYAR
  498.       IF (KDESIGN.EQ.1) THEN
  499. *      IF (NP.NE.3) RETURN
  500.       S=0.d0
  501.       do 120 I=1,NP
  502.       S=S+abs(UMI(KP,1)*HR(K,I,1)+UMI(KP,2)*HR(K,I,2))/UM(KP)
  503.   120 CONTINUE
  504.       HMK=2.d0/S
  505.       S=0.d0
  506.       do 121 I=1,NP
  507.       S=S+abs(UPI(KP,1)*HR(K,I,1)+UPI(KP,2)*HR(K,I,2))/UP(KP)
  508.   121 CONTINUE
  509.       HPK=2.d0/S
  510. C DESIGN ORIGINAL
  511.       ELSEIF (KDESIGN.EQ.0) THEN
  512.       HMK=2.D0*UM(KP)/BM(KP)
  513.       HPK=2.D0*UP(KP)/BP(KP)
  514.  
  515.       ELSEIF (KDESIGN.EQ.2) THEN
  516.       do 777 i=1,NP
  517.       x(i)=XCOOR((LE(i,K)-1)*(IES+1)+1)
  518.       y(i)=XCOOR((LE(i,K)-1)*(IES+1)+2)
  519.   777 CONTINUE
  520.  
  521.       DJ=(x(2)-x(1)+(x(3)-x(4))*(NP-3))*(y(3)-y(1)+(y(4)-y(2))*(NP-3))
  522.      *  -(y(2)-y(1)+(y(3)-y(4))*(NP-3))*(x(3)-x(1)+(x(4)-x(2))*(NP-3))
  523.  
  524.       IF (NP.EQ.3) lll=1.d0
  525.       IF (NP.EQ.4) lll=4.d0
  526.       DJ=DJ/(lll**2)
  527.  
  528.       b(1)=(UMI(KP,1)*(y(3)-y(1)+(y(4)-y(2))*(NP-3))-UMI(KP,2)*(x(3)-
  529.      * x(1)+(x(4)-x(2))*(NP-3)))/DJ/lll
  530.       b(2)=(-UMI(KP,1)*(y(2)-y(1)+(y(3)-y(4))*(NP-3))+UMI(KP,2)*(x(2)-
  531.      * x(1)+(x(3)-x(4))*(NP-3)))/DJ/lll
  532.  
  533.       p=1
  534.       s=0.d0
  535.       do 306 kk=1,2
  536.       s=s+(abs(b(kk)))**p
  537.   306 CONTINUE
  538.       Nm=s**(1.d0/p)
  539.       HMK=2.d0*UM(KP)/Nm
  540. C Calcul de HPK:
  541.  
  542.       b(1)=(UPI(KP,1)*(y(3)-y(1)+(y(4)-y(2))*(NP-3))-UPI(KP,2)*(x(3)-
  543.      * x(1)+(x(4)-x(2))*(NP-3)))/DJ/lll
  544.       b(2)=(-UPI(KP,1)*(y(2)-y(1)+(y(3)-y(4))*(NP-3))+UPI(KP,2)*(x(2)-
  545.      * x(1)+(x(3)-x(4))*(NP-3)))/DJ/lll
  546.  
  547.       s=0.d0
  548.       do 308 kk=1,2
  549.       s=s+(abs(b(kk)))**p
  550.   308 CONTINUE
  551.       Nm=s**(1.d0/p)
  552.       HPK=2.d0*UP(KP)/Nm
  553.  
  554.       ELSEIF (KDESIGN.EQ.3) THEN
  555.       do 778 i=1,NP
  556.       x(i)=XCOOR((LE(i,K)-1)*(IES+1)+1)
  557.       y(i)=XCOOR((LE(i,K)-1)*(IES+1)+2)
  558.   778 CONTINUE
  559.  
  560.       DJ=(x(2)-x(1)+(x(3)-x(4))*(NP-3))*(y(3)-y(1)+(y(4)-y(2))*(NP-3))
  561.      *  -(y(2)-y(1)+(y(3)-y(4))*(NP-3))*(x(3)-x(1)+(x(4)-x(2))*(NP-3))
  562.  
  563.       IF (NP.EQ.3) lll=1.d0
  564.       IF (NP.EQ.4) lll=4.d0
  565.       DJ=DJ/(lll**2)
  566.  
  567.       b(1)=(UMI(KP,1)*(y(3)-y(1)+(y(4)-y(2))*(NP-3))-UMI(KP,2)*(x(3)-
  568.      * x(1)+(x(4)-x(2))*(NP-3)))/DJ/lll
  569.       b(2)=(-UMI(KP,1)*(y(2)-y(1)+(y(3)-y(4))*(NP-3))+UMI(KP,2)*(x(2)-
  570.      * x(1)+(x(3)-x(4))*(NP-3)))/DJ/lll
  571.  
  572.       p=2
  573.       s=0.d0
  574.       do 309 kk=1,2
  575.       s=s+(abs(b(kk)))**p
  576.   309 CONTINUE
  577.       Nm=s**(1.d0/p)
  578.       HMK=2.d0*UM(KP)/Nm
  579. C Calcul de HPK:
  580.  
  581.       b(1)=(UPI(KP,1)*(y(3)-y(1)+(y(4)-y(2))*(NP-3))-UPI(KP,2)*(x(3)-
  582.      * x(1)+(x(4)-x(2))*(NP-3)))/DJ/lll
  583.       b(2)=(-UPI(KP,1)*(y(2)-y(1)+(y(3)-y(4))*(NP-3))+UPI(KP,2)*(x(2)-
  584.      * x(1)+(x(3)-x(4))*(NP-3)))/DJ/lll
  585.  
  586.       s=0.d0
  587.       do 310 kk=1,2
  588.       s=s+(abs(b(kk)))**p
  589.   310 CONTINUE
  590.       Nm=s**(1.d0/p)
  591.       HPK=2.d0*UP(KP)/Nm
  592.  
  593.       ELSEIF (KDESIGN.EQ.4) THEN
  594.       do 779 i=1,NP
  595.       x(i)=XCOOR((LE(i,K)-1)*(IES+1)+1)
  596.       y(i)=XCOOR((LE(i,K)-1)*(IES+1)+2)
  597.   779 CONTINUE
  598.  
  599. * ZZ=1 :HMK=HPK
  600. * ZZ=2 :HMK<>HPK
  601.  
  602.       do 321 i=1,NP
  603.       kkk(i)=Aire(KP)*(UMI(KP,1)*HR(K,i,1)+UMI(KP,2)*HR(K,i,2))
  604.   321 CONTINUE
  605.  
  606.       s0=0.d0
  607.       s1=0.d0
  608.       do 322 i=1,NP
  609.       s0=s0+min(0.d0,kkk(i))
  610.       s1=s1+max(0.d0,kkk(i))
  611.   322 CONTINUE
  612.  
  613.       s2=0.d0
  614.       s3=0.d0
  615.       do 323 i=1,NP
  616.       s2=s2-min(0.d0,kkk(i))*x(i)/s0+max(0.d0,kkk(i))*x(i)/s1
  617.       s3=s3-min(0.d0,kkk(i))*y(i)/s0+max(0.d0,kkk(i))*y(i)/s1
  618.   323 CONTINUE
  619.       HMK=sqrt(s2**2+s3**2)
  620.       HPK=HMK
  621.  
  622.       IF (zz.eq.2) THEN
  623.       do 324 i=1,NP
  624.       kkk(i)=Aire(KP)*(UPI(KP,1)*HR(K,i,1)+UPI(KP,2)*HR(K,i,2))
  625.   324 CONTINUE
  626.  
  627.       s0=0.d0
  628.       s1=0.d0
  629.       do 325 i=1,NP
  630.       s0=s0+min(0.d0,kkk(i))
  631.       s1=s1+max(0.d0,kkk(i))
  632.   325 CONTINUE
  633.  
  634.       s2=0.d0
  635.       s3=0.d0
  636.       do 3266 i=1,NP
  637.       s2=s2-min(0.d0,kkk(i))*x(i)/s0+max(0.d0,kkk(i))*x(i)/s1
  638.       s3=s3-min(0.d0,kkk(i))*y(i)/s0+max(0.d0,kkk(i))*y(i)/s1
  639.  3266 CONTINUE
  640.       HPK=sqrt(s2**2+s3**2)
  641.       ENDIF
  642.  
  643.       ELSEIF (KDESIGN.EQ.5) THEN
  644.       do 802 i=1,NP
  645.       n(1,i)=abs(2.d0*Aire(KP)*HR(K,i,1))
  646.       n(2,i)=abs(2.d0*Aire(KP)*HR(K,i,2))
  647.       MMAX(i)=n(1,i)
  648.       IF (n(1,i).lt.n(2,i)) MMAX(i)=n(2,i)
  649.   802 CONTINUE
  650.  
  651.       HMK=0.d0
  652.       do 803 i=1,NP
  653.       IF (MMAX(i).gt.HMK) HMK=MMAX(i)
  654.   803 CONTINUE
  655.  
  656.       ENDIF
  657.  
  658.       IF (HMK.LT.HMIN) HMIN=HMK
  659.       IF (HMK.GT.HMAX) HMAX=HMK
  660.  
  661.       XMB(KP)=HMK
  662.       ALFA=UM(KP)*HMK/ALFT(KP)/2.D0
  663.       AKSI=ALFA/3.D0
  664.       IF(ALFA.GT.3.D0)AKSI=1.D0
  665.       CCT=AKSI*HMK/2.d0/UM(KP)
  666.  
  667. * si on a le meme h
  668.       HPK=HMK
  669.       ALFA=UM(KP)*HPK/ALFT(KP)/2.D0
  670.       AKSI=ALFA/3.D0
  671.       IF(ALFA.GT.3.D0)AKSI=1.D0
  672.       CCP=AKSI*HPK/2.d0/UP(KP)
  673.       CPT=CCP-CCT
  674.       CC2=0.D0
  675.       IF(CPT.GE.0.D0)CC2=CPT
  676.  
  677.       CXT(KP)=UMI(KP,1)*CCT
  678.       CYT(KP)=UMI(KP,2)*CCT
  679.       CXY(KP)=CCT
  680.  
  681.       DXT(KP)=UPI(KP,1)*UPI(KP,1)*CC2
  682.       DYT(KP)=UPI(KP,2)*UPI(KP,2)*CC2
  683.       DXY(KP)=UPI(KP,1)*UPI(KP,2)*CC2
  684.  
  685.  7108 CONTINUE
  686.  
  687.       ELSEIF(IDCENN.EQ.3)THEN
  688.       DO 71018 K=KDEB,KFIN
  689.       KP=K-KDEB+1
  690. C DESIGN TEZDUYAR
  691.       IF (KDESIGN.EQ.1) THEN
  692. *      IF (NP.NE.3) RETURN
  693.       S=0.d0
  694.       do 122 I=1,NP
  695.       S=S+abs(UMI(KP,1)*HR(K,I,1)+UMI(KP,2)*HR(K,I,2))/UM(KP)
  696.   122 CONTINUE
  697.       HMK=2.d0/S
  698. C DESIGN ORIGINAL
  699.       ELSEIF (KDESIGN.EQ.0) THEN
  700.       HMK=2.D0*UM(KP)/BM(KP)
  701.  
  702.       ELSEIF (KDESIGN.EQ.2) THEN
  703.       do 780 i=1,NP
  704.       x(i)=XCOOR((LE(i,K)-1)*(IES+1)+1)
  705.       y(i)=XCOOR((LE(i,K)-1)*(IES+1)+2)
  706.   780 CONTINUE
  707.  
  708.       DJ=(x(2)-x(1)+(x(3)-x(4))*(NP-3))*(y(3)-y(1)+(y(4)-y(2))*(NP-3))
  709.      *  -(y(2)-y(1)+(y(3)-y(4))*(NP-3))*(x(3)-x(1)+(x(4)-x(2))*(NP-3))
  710.  
  711.       IF (NP.EQ.3) lll=1.d0
  712.       IF (NP.EQ.4) lll=4.d0
  713.       DJ=DJ/(lll**2)
  714.  
  715.       b(1)=(UMI(KP,1)*(y(3)-y(1)+(y(4)-y(2))*(NP-3))-UMI(KP,2)*(x(3)-
  716.      * x(1)+(x(4)-x(2))*(NP-3)))/DJ/lll
  717.       b(2)=(-UMI(KP,1)*(y(2)-y(1)+(y(3)-y(4))*(NP-3))+UMI(KP,2)*(x(2)-
  718.      * x(1)+(x(3)-x(4))*(NP-3)))/DJ/lll
  719.  
  720.       p=1
  721.       s=0.d0
  722.       do 311 kk=1,2
  723.       s=s+(abs(b(kk)))**p
  724.   311 CONTINUE
  725.       Nm=s**(1.d0/p)
  726.       HMK=2.d0*UM(KP)/Nm
  727.  
  728.  
  729.       ELSEIF (KDESIGN.EQ.3) THEN
  730.       do 781 i=1,NP
  731.       x(i)=XCOOR((LE(i,K)-1)*(IES+1)+1)
  732.       y(i)=XCOOR((LE(i,K)-1)*(IES+1)+2)
  733.   781 CONTINUE
  734.  
  735.       DJ=(x(2)-x(1)+(x(3)-x(4))*(NP-3))*(y(3)-y(1)+(y(4)-y(2))*(NP-3))
  736.      *  -(y(2)-y(1)+(y(3)-y(4))*(NP-3))*(x(3)-x(1)+(x(4)-x(2))*(NP-3))
  737.  
  738.       IF (NP.EQ.3) lll=1.d0
  739.       IF (NP.EQ.4) lll=4.d0
  740.       DJ=DJ/(lll**2)
  741.  
  742.       b(1)=(UMI(KP,1)*(y(3)-y(1)+(y(4)-y(2))*(NP-3))-UMI(KP,2)*(x(3)-
  743.      * x(1)+(x(4)-x(2))*(NP-3)))/DJ/lll
  744.       b(2)=(-UMI(KP,1)*(y(2)-y(1)+(y(3)-y(4))*(NP-3))+UMI(KP,2)*(x(2)-
  745.      * x(1)+(x(3)-x(4))*(NP-3)))/DJ/lll
  746.  
  747.       p=2
  748.       s=0.d0
  749.       do 312 kk=1,2
  750.       s=s+(abs(b(kk)))**p
  751.   312 CONTINUE
  752.       Nm=s**(1.d0/p)
  753.       HMK=2.d0*UM(KP)/Nm
  754.  
  755.       ELSEIF (KDESIGN.EQ.4) THEN
  756.       do 782 i=1,NP
  757.       x(i)=XCOOR((LE(i,K)-1)*(IES+1)+1)
  758.       y(i)=XCOOR((LE(i,K)-1)*(IES+1)+2)
  759.   782 CONTINUE
  760.  
  761.       do 326 i=1,NP
  762.       kkk(i)=Aire(KP)*(UMI(KP,1)*HR(K,i,1)+UMI(KP,2)*HR(K,i,2))
  763.   326 CONTINUE
  764.  
  765.       s0=0.d0
  766.       s1=0.d0
  767.       do 3277 i=1,KP
  768.       s0=s0+min(0.d0,kkk(i))
  769.       s1=s1+max(0.d0,kkk(i))
  770.  3277 CONTINUE
  771.  
  772.       s2=0.d0
  773.       s3=0.d0
  774.       do 3288 i=1,KP
  775.       s2=s2-min(0.d0,kkk(i))*x(i)/s0+max(0.d0,kkk(i))*x(i)/s1
  776.       s3=s3-min(0.d0,kkk(i))*y(i)/s0+max(0.d0,kkk(i))*y(i)/s1
  777.  3288 CONTINUE
  778.       HMK=sqrt(s2**2+s3**2)
  779.  
  780.       ELSEIF (KDESIGN.EQ.5) THEN
  781.       do 804 i=1,NP
  782.       n(1,i)=abs(2.d0*Aire(KP)*HR(K,i,1))
  783.       n(2,i)=abs(2.d0*Aire(KP)*HR(K,i,2))
  784.       MMAX(i)=n(1,i)
  785.       IF (n(1,i).lt.n(2,i)) MMAX(i)=n(2,i)
  786.   804 CONTINUE
  787.  
  788.       HMK=0.d0
  789.       do 805 i=1,NP
  790.       IF (MMAX(i).gt.HMK) HMK=MMAX(i)
  791.   805 CONTINUE
  792.  
  793.       ENDIF
  794.  
  795.       IF (HMK.LT.HMIN) HMIN=HMK
  796.       IF (HMK.GT.HMAX) HMAX=HMK
  797.  
  798.       XMB(KP)=HMK
  799.       ALFA=UM(KP)*HMK/ALFT(KP)/2.D0
  800.       AKSI=ALFA/3.D0
  801.       IF(ALFA.GT.3.D0)AKSI=1.D0
  802.       CCT=AKSI*HMK/2.d0/UM(KP)
  803.  
  804.       CXT(KP)=UMI(KP,1)*CCT
  805.       CYT(KP)=UMI(KP,2)*CCT
  806.       CXY(KP)=CCT
  807.  
  808.       DXT(KP)=0.D0
  809.       DYT(KP)=0.D0
  810.       DXY(KP)=0.D0
  811.  
  812. 71018 CONTINUE
  813.  
  814.       ELSEIF(IDCENN.EQ.4)THEN
  815.       DT19=DTM1*0.5D0
  816.       DO 71009 K=KDEB,KFIN
  817.       KP=K-KDEB+1
  818.       XMB(KP)=XMH(KP)
  819.       CXT(KP)=UMI(KP,1)*DT19
  820.       CYT(KP)=UMI(KP,2)*DT19
  821.       CXY(KP)=DT19
  822.  
  823.       DXT(KP)=0.D0
  824.       DYT(KP)=0.D0
  825.       DXY(KP)=0.D0
  826.  
  827. 71009 CONTINUE
  828.  
  829.       ENDIF
  830.  
  831.       ENDIF
  832. C           ***********************
  833. C
  834. C     AVANT CALCUL DT
  835. C
  836.       IF(IKOMP.EQ.0)THEN
  837.  
  838. C*IBMDIR* PREFER SCALAR
  839.       DO 7010 K=KDEB,KFIN
  840.       KP=K-KDEB+1
  841.       DT0=DT
  842.  
  843.       DT1=2.D0/
  844.      & ( (UMI(KP,1)*UMI(KP,1))/(ALFT(KP)+CXT(KP))
  845.      & + (UMI(KP,2)*UMI(KP,2))/(ALFT(KP)+CYT(KP)) )
  846.  
  847.       DT2=0.5D0/
  848.      & ( (ALFT(KP)+CXT(KP))*AL2(KP)
  849.      &  +(ALFT(KP)+CYT(KP))*AH2(KP) )
  850.  
  851.  
  852.       IF(DT1.LT.DT)DT=DT1
  853.       IF(DT2.LT.DT)DT=DT2
  854. C     IF(DT3.LT.DT)DT=DT3
  855.       IF(DT.NE.DT0) THEN
  856.          DTT1=DT1
  857.          DTT2=DT2
  858.          DTT3=0.D0
  859.          DIAEL=XMB(KP)
  860.          NUEL=K
  861.       END IF
  862.       CXY(KP)=CXY(KP)*CD
  863.  7010 CONTINUE
  864.  
  865.       ELSEIF(IKOMP.EQ.1)THEN
  866.  
  867. C*IBMDIR* PREFER SCALAR
  868.       DO 7110 K=KDEB,KFIN
  869.       KP=K-KDEB+1
  870.       DT0=DT
  871.  
  872.       DT1=2.D0/
  873.      & ( (UMI(KP,1)*UMI(KP,1))/(ALFT(KP)+CXT(KP)*UMI(KP,1)+DXT(KP))
  874.      & + (UMI(KP,2)*UMI(KP,2))/(ALFT(KP)+CYT(KP)*UMI(KP,2)+DYT(KP)) )
  875.  
  876.       DT2=0.5D0/
  877.      & ( (ALFT(KP)+CXT(KP)*UMI(KP,1)+DXT(KP))*AL2(KP)
  878.      &  +(ALFT(KP)+CYT(KP)*UMI(KP,2)+DYT(KP))*AH2(KP) )
  879.  
  880.       IF(DT1.LT.DT)DT=DT1
  881.       IF(DT2.LT.DT)DT=DT2
  882. C     IF(DT3.LT.DT)DT=DT3
  883.       IF(DT.NE.DT0) THEN
  884.          DTT1=DT1
  885.          DTT2=DT2
  886.          DTT3=0.D0
  887.          DIAEL=XMB(KP)
  888.          NUEL=K
  889.       END IF
  890.       CXY(KP)=CXY(KP)*CD
  891.  7110 CONTINUE
  892.  
  893.       ENDIF
  894. C Le coeur du calcul ...
  895.  
  896.       IF(IKOMP.EQ.0)THEN
  897.  
  898.       DO 81086 I=1,NP
  899.       DO 81085 J= 1,NP
  900. C*IBMDIR* PREFER VECTOR
  901.  
  902.       DO 7014 K=KDEB,KFIN
  903.       KP=K-KDEB+1
  904.  
  905.       ZVGG=AIRE(KP)*(
  906.      &     HR(K,I,1)*HR(K,J,1)*CXT(KP)
  907.      &+    HR(K,I,1)*HR(K,J,2)*CXY(KP)
  908.      &+    HR(K,I,2)*HR(K,J,1)*CXY(KP)
  909.      &+    HR(K,I,2)*HR(K,J,2)*CYT(KP)
  910.      &+    (CXT(KP)*AL2(KP)+CYT(KP)*AH2(KP))/12.D0
  911.      &     *VGGT(J,I)*QUA4  )
  912.  
  913. C    &+    HR(K,I,2)*HR(K,J,2)*CYT(KP) )
  914. C    &+    (CXT(KP)+CYT(KP))/24.D0*VGGT(J,I)*QUA4
  915.  
  916.       ZVGT=AIRE(KP)*(
  917.      &     HR(K,I,1)*HR(K,J,1)*ALFT(KP)
  918.      &+    HR(K,I,2)*HR(K,J,2)*ALFT(KP)
  919.      &+    (ALFT(KP)*AL2(KP)+ALFT(KP)*AH2(KP))/12.D0
  920.      &     *VGGT(J,I)*QUA4  )
  921.  
  922. C    &+    HR(K,I,2)*HR(K,J,2)*ALFT(KP) )
  923. C    &+    ALFT(KP)/12.D0*VGGT(J,I)*QUA4
  924.  
  925. C?    V2=(UIX(KP,I)*HR(K,J,1)+UIY(KP,I)*HR(K,J,2))*DR(KP,I)
  926.       V2=UMI(KP,1)*DR(KP,I)*HR(K,J,1)+UMI(KP,2)*DR(KP,I)*HR(K,J,2)
  927.  
  928.       SBF(KP,I)=SBF(KP,I)+TETAC(KP,J)*(ZVGG+V2)+ TETAD(KP,J)*ZVGT
  929.  
  930. 7014  CONTINUE
  931. 81085 CONTINUE
  932. 81086 CONTINUE
  933.       ELSEIF(IKOMP.EQ.1)THEN
  934.  
  935.       DO 81088 I=1,NP
  936.       DO 81087 J= 1,NP
  937. C*IBMDIR* PREFER VECTOR
  938.       DO 7015 K=KDEB,KFIN
  939.       KP=K-KDEB+1
  940.       ZVGG=AIRE(KP)*(
  941.      &     HR(K,I,1)*HR(K,J,1)*(CXT(KP)*UIX(KP,J)+DXT(KP))
  942.      &+    HR(K,I,1)*HR(K,J,2)
  943.      &     *(CXY(KP)*UMI(KP,1)*UIY(KP,J)+DXY(KP))
  944.      &+    HR(K,I,2)*HR(K,J,1)
  945.      &     *(CXY(KP)*UMI(KP,2)*UIX(KP,J)+DXY(KP))
  946.      &+    HR(K,I,2)*HR(K,J,2)*(CYT(KP)*UIY(KP,J)+DYT(KP))
  947.      &+    ((CXT(KP)*UIX(KP,J)+DXT(KP))*AL2(KP)
  948.      &      +(CYT(KP)*UIY(KP,J)+DYT(KP))*AH2(KP))/12.D0
  949.      &     *VGGT(J,I)*QUA4  )
  950.  
  951. C    &+    HR(K,I,2)*HR(K,J,2)*CYT(KP) )
  952. C    &+    (CXT(KP)+CYT(KP))/24.D0*VGGT(J,I)*QUA4
  953.  
  954.       ZVGT=AIRE(KP)*(
  955.      &     HR(K,I,1)*HR(K,J,1)*ALFT(KP)
  956.      &+    HR(K,I,2)*HR(K,J,2)*ALFT(KP)
  957.      &+    (ALFT(KP)*AL2(KP)+ALFT(KP)*AH2(KP))/12.D0
  958.      &     *VGGT(J,I)*QUA4  )
  959.  
  960. C    &+    HR(K,I,2)*HR(K,J,2)*ALFT(KP) )
  961. C    &+    ALFT(KP)/12.D0*VGGT(J,I)*QUA4
  962.  
  963.       V2=(UIX(KP,J)*HR(K,J,1)+UIY(KP,J)*HR(K,J,2))*DR(KP,I)
  964.  
  965.       SBF(KP,I)=SBF(KP,I)+TETAC(KP,J)*(ZVGG+V2)+ TETAD(KP,J)*ZVGT
  966.  
  967.  7015 CONTINUE
  968. 81087 CONTINUE
  969. 81088 CONTINUE
  970.  
  971.       ENDIF
  972. C
  973. C     Fin de l'accumulation dans SBF.
  974. C     On ajoute ces incréments G.
  975. C
  976.       DO 81089 I=1,NP
  977. C*IBMDIR* PREFER VECTOR
  978.       DO 7017 K=KDEB,KFIN
  979.       KP=K-KDEB+1
  980.       NF=IPADS(LE(I,K))
  981.       G(NF) = G(NF)-SBF(KP,I)
  982.  7017 CONTINUE
  983. 81089 CONTINUE
  984.  
  985.  7001 CONTINUE
  986.  
  987. C     WRITE(6,*)' G DANS YCTSCL '
  988. C     WRITE(6,1984)(M,G(M),M=1,NPTD)
  989.  1984 FORMAT(7(1X,I4,2X,1PE11.4))
  990.  
  991. C     CALL ARRET(0)
  992.       IPAS=1
  993. C
  994. C
  995. C      PRINT *,'HMIN = ', HMIN, 'HMAX = ', HMAX
  996. C
  997.       RETURN
  998.  
  999. C                           ********
  1000. C                           *  3D  *
  1001. C                           ********
  1002.  
  1003.  10   CONTINUE
  1004.  
  1005.       IF(IPAS.EQ.0)CALL CALHRH(QGGT,Q1,Q2,Q3,IES)
  1006.       CUB8=0.D0
  1007.       IF(NP.EQ.8)CUB8=1.D0
  1008.  
  1009. C
  1010. C Calcul du nombre de paquets de LRV éléments
  1011. C
  1012.       NNN=MOD(NEL,LRV)
  1013.       IF(NNN.EQ.0) NPACK=NEL/LRV
  1014.       IF(NNN.NE.0) NPACK=1+(NEL-NNN)/LRV
  1015.       KPACKD=1
  1016.       KPACKF=NPACK
  1017. C
  1018. C ******* BOUCLE SUR LES PAQUETS DE LRV ELEMENTS **********
  1019. C
  1020.       DO 8001 KPACK=KPACKD,KPACKF
  1021. C
  1022. C ======= A L'INTERIEUR DE CHAQUE PAQUET DE LRV ELEMENTS =======
  1023. C
  1024. C 1. Calcul des limites du paquet courant.
  1025.       KDEB=1+(KPACK-1)*LRV
  1026.       KFIN=MIN(NEL,KDEB+LRV-1)
  1027. C
  1028.       DO 8002 K=KDEB,KFIN
  1029.       KP=K-KDEB+1
  1030.       NK=K+K0
  1031.       NK1=(1-IND1)*(NK-1)+1
  1032.       ALFT(KP)=ALFE(NK1)+XPETIT
  1033.       AIRE(KP)=VOLU(NK)
  1034.       AL(KP)=COTE(NK,1)
  1035.       AH(KP)=COTE(NK,2)
  1036.       AP(KP)=COTE(NK,3)
  1037.       CFM(KP)=AL(KP)*AH(KP)/AP(KP)+AL(KP)*AP(KP)/AH(KP)+
  1038.      &        AP(KP)*AH(KP)/AL(KP)
  1039. C     CF1(KP)=AL(KP)*AH(KP)/AP(KP)
  1040. C     CF2(KP)=AL(KP)*AP(KP)/AH(KP)
  1041. C     CF3(KP)=AP(KP)*AH(KP)/AL(KP)
  1042.       XMH(KP)=(AL(KP)+AH(KP)+AP(KP))/3.D0
  1043.       AL2(KP)=1.D0/AL(KP)/AL(KP)
  1044.       AH2(KP)=1.D0/AH(KP)/AH(KP)
  1045.       AP2(KP)=1.D0/AP(KP)/AP(KP)
  1046.  8002 CONTINUE
  1047.       IF((IKOMP.EQ.0.AND.IKAS.EQ.5).OR.
  1048.      &(IKOMP.EQ.1.AND.IKAS.EQ.6))THEN
  1049.       DO 8003 K=KDEB,KFIN
  1050.       KP=K-KDEB+1
  1051.       NK=K+K0
  1052.       ALFT(KP)=ALFT(KP)+ALT(NK)/SGT
  1053.  8003 CONTINUE
  1054.       ENDIF
  1055.  
  1056. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  1057.  
  1058. C Initialisation des UMI avant accumulation
  1059.       DO 8005 K=KDEB,KFIN
  1060.       KP=K-KDEB+1
  1061.       UMI(KP,1)=XPETIT
  1062.       UMI(KP,2)=XPETIT
  1063.       UMI(KP,3)=XPETIT
  1064.       UPI(KP,1)=XPETIT
  1065.       UPI(KP,2)=XPETIT
  1066.       UPI(KP,3)=XPETIT
  1067.       GRADT(KP,1)=XPETIT
  1068.       GRADT(KP,2)=XPETIT
  1069.       GRADT(KP,3)=XPETIT
  1070.  8005 CONTINUE
  1071.  
  1072.       DO 81090 I=1,NP
  1073. C*IBMDIR* PREFER VECTOR
  1074.       DO 8006 K=KDEB,KFIN
  1075.       KP=K-KDEB+1
  1076.       NF=IPADI(LE(I,K))
  1077.       NFU=IPADS(LE(I,K))
  1078.       NFU=(1-INDU)*(NFU-1)+1
  1079.       DR(KP,I)=DRR(I,K)
  1080.       UIX(KP,I)=UN(NFU,1)
  1081.       UIY(KP,I)=UN(NFU,2)
  1082.       UIZ(KP,I)=UN(NFU,3)
  1083.       UMI(KP,1)=UMI(KP,1)+UN(NFU,1)*DR(KP,I)
  1084.       UMI(KP,2)=UMI(KP,2)+UN(NFU,2)*DR(KP,I)
  1085.       UMI(KP,3)=UMI(KP,3)+UN(NFU,3)*DR(KP,I)
  1086.  
  1087.       TETAC(KP,I)=HRN(NF)
  1088.       TETAD(KP,I)=TN(NF)
  1089.       GRADT(KP,1)=GRADT(KP,1)+HR(K,I,1)*TETAC(KP,I)
  1090.       GRADT(KP,2)=GRADT(KP,2)+HR(K,I,2)*TETAC(KP,I)
  1091.       GRADT(KP,3)=GRADT(KP,3)+HR(K,I,3)*TETAC(KP,I)
  1092.  8006 CONTINUE
  1093. 81090 CONTINUE
  1094.  
  1095. C
  1096. C Initialisation de la variable d'accumulation SBF au terme source
  1097. C                                                                      M
  1098. C       write(6,*)' IKomp,ikas=',IKomp,ikas
  1099. C       write(6,*)' IKS,IND1,INDU=',IKS,IND1,INDU
  1100.  
  1101.       IF(IKOMP.EQ.0)THEN
  1102.  
  1103.          DO 80021 K=KDEB,KFIN
  1104.          KP=K-KDEB+1
  1105.          NK=K+K0
  1106.          NKS=(1-IKS)*(NK-1)+1
  1107.          QT(KP)=QE(NKS)
  1108. 80021    CONTINUE
  1109.  
  1110.          DO 81091 I=1,NP
  1111. C*IBMDIR* PREFER VECTOR
  1112.          DO 80062 K=KDEB,KFIN
  1113.          KP=K-KDEB+1
  1114.          SBF(KP,I)=-QT(KP)*DR(KP,I)
  1115. 80062    CONTINUE
  1116. 81091    CONTINUE
  1117.  
  1118.       ELSEIF(IKOMP.EQ.1)THEN
  1119.  
  1120.          DO 80023 K=KDEB,KFIN
  1121.          KP=K-KDEB+1
  1122.          NK=K+K0
  1123.          NKS=(1-IKS)*(NK-1)+1
  1124.          QT(KP)=QE(NKS)
  1125. 80023    CONTINUE
  1126.  
  1127.          DO 81092 I=1,NP
  1128. C*IBMDIR* PREFER VECTOR
  1129.          DO 80064 K=KDEB,KFIN
  1130.          KP=K-KDEB+1
  1131.          SBF(KP,I)=-QT(KP)*DR(KP,I)
  1132. 80064    CONTINUE
  1133. 81092    CONTINUE
  1134.  
  1135.       ENDIF
  1136.  
  1137.       DO 8007 K=KDEB,KFIN
  1138.       KP=K-KDEB+1
  1139.       UMI(KP,1)=UMI(KP,1)/AIRE(KP)
  1140.       UMI(KP,2)=UMI(KP,2)/AIRE(KP)
  1141.       UMI(KP,3)=UMI(KP,3)/AIRE(KP)
  1142.       UM(KP)=UMI(KP,1)*UMI(KP,1)+UMI(KP,2)*UMI(KP,2)
  1143.      &      +UMI(KP,3)*UMI(KP,3)
  1144.       UM(KP)=SQRT(UM(KP))
  1145.       A=GRADT(KP,1)*GRADT(KP,1)+GRADT(KP,2)*GRADT(KP,2)
  1146.      & +GRADT(KP,3)*GRADT(KP,3)
  1147.       A=SQRT(A)+XPETIT
  1148.       GRADT(KP,1)=GRADT(KP,1)/A
  1149.       GRADT(KP,2)=GRADT(KP,2)/A
  1150.       GRADT(KP,3)=GRADT(KP,3)/A
  1151.       UP(KP)=GRADT(KP,1)*UMI(KP,1)+GRADT(KP,2)*UMI(KP,2)
  1152.      &      +GRADT(KP,3)*UMI(KP,3)
  1153.       UPI(KP,1)=UP(KP)*GRADT(KP,1)
  1154.       UPI(KP,2)=UP(KP)*GRADT(KP,2)
  1155.       UPI(KP,3)=UP(KP)*GRADT(KP,3)
  1156.  8007 CONTINUE
  1157.  
  1158. *      DO 80074 K=KDEB,KFIN
  1159. *      KP=K-KDEB+1
  1160. *      BMX=UMI(KP,1)/AL(KP)
  1161. *      BMY=UMI(KP,2)/AH(KP)
  1162. *      BMZ=UMI(KP,3)/AP(KP)
  1163. *      BM(KP)=BMX*BMX+BMY*BMY+BMZ*BMZ
  1164. *      BM(KP)=SQRT(BM(KP))+XPETIT
  1165. *      BPX=UPI(KP,1)/AL(KP)
  1166. *      BPY=UPI(KP,2)/AH(KP)
  1167. *      BPZ=UPI(KP,3)/AP(KP)
  1168. *      BP(KP)=BPX*BPX+BPY*BPY+BPZ*BPZ
  1169. *      BP(KP)=SQRT(BP(KP))+XPETIT
  1170. *80074 CONTINUE
  1171.  
  1172. C
  1173. C     AVANT DECENTREMENT
  1174. C
  1175.  
  1176.       IF(IKOMP.EQ.0)THEN
  1177.  
  1178. C ALFA,AKSI,CCU,CCT utilisées seulement dans la boucle 8008
  1179.       IF(IDCENN.EQ.1)THEN
  1180.       DO 80081 K=KDEB,KFIN
  1181.       KP=K-KDEB+1
  1182.       XMB(KP)=XMH(KP)
  1183.       CXT(KP)=0.D0
  1184.       CYT(KP)=0.D0
  1185.       CZT(KP)=0.D0
  1186.       CXY(KP)=0.D0
  1187.       CXZ(KP)=0.D0
  1188.       CYZ(KP)=0.D0
  1189. 80081 CONTINUE
  1190.  
  1191.       ELSEIF(IDCENN.EQ.2)THEN
  1192.       DO 8008 K=KDEB,KFIN
  1193.       KP=K-KDEB+1
  1194.       HMK=2.D0*UM(KP)/BM(KP)
  1195.       XMB(KP)=HMK
  1196.       ALFA=UM(KP)*HMK/ALFT(KP)/2.D0
  1197.       AKSI=ALFA/3.D0
  1198.       IF(ALFA.GT.3.D0)AKSI=1.D0
  1199.       CCT=AKSI/BM(KP)
  1200.  
  1201.       HPK=2.D0*UP(KP)/BP(KP)
  1202.       ALFA=UP(KP)*HPK/ALFT(KP)/2.D0
  1203.       AKSI=ALFA/3.D0
  1204.       IF(ALFA.GT.3.D0)AKSI=1.D0
  1205.       CCP=AKSI/BP(KP)
  1206.       CPT=CCP-CCT
  1207.       CC2=0.D0
  1208.       IF(CPT.GE.0.D0)CC2=CPT
  1209.  
  1210.       CXT(KP)=(UMI(KP,1)*UMI(KP,1)*CCT+UPI(KP,1)*UPI(KP,1)*CC2)
  1211.       CYT(KP)=(UMI(KP,2)*UMI(KP,2)*CCT+UPI(KP,2)*UPI(KP,2)*CC2)
  1212.       CZT(KP)=(UMI(KP,3)*UMI(KP,3)*CCT+UPI(KP,3)*UPI(KP,3)*CC2)
  1213.       CXY(KP)=(UMI(KP,1)*UMI(KP,2)*CCT+UPI(KP,1)*UPI(KP,2)*CC2)
  1214.       CXZ(KP)=(UMI(KP,1)*UMI(KP,3)*CCT+UPI(KP,1)*UPI(KP,3)*CC2)
  1215.       CYZ(KP)=(UMI(KP,2)*UMI(KP,3)*CCT+UPI(KP,2)*UPI(KP,3)*CC2)
  1216.  
  1217.  8008 CONTINUE
  1218.  
  1219.       ELSEIF(IDCENN.EQ.3)THEN
  1220.       DO 8018 K=KDEB,KFIN
  1221.       KP=K-KDEB+1
  1222.       HMK=2.D0*UM(KP)/BM(KP)
  1223.       XMB(KP)=HMK
  1224.       ALFA=UM(KP)*HMK/ALFT(KP)/2.D0
  1225.       AKSI=ALFA/3.D0
  1226.       IF(ALFA.GT.3.D0)AKSI=1.D0
  1227.       CCT=AKSI/BM(KP)
  1228.  
  1229.       CXT(KP)=UMI(KP,1)*UMI(KP,1)*CCT
  1230.       CYT(KP)=UMI(KP,2)*UMI(KP,2)*CCT
  1231.       CZT(KP)=UMI(KP,3)*UMI(KP,3)*CCT
  1232.       CXY(KP)=UMI(KP,1)*UMI(KP,2)*CCT
  1233.       CXZ(KP)=UMI(KP,1)*UMI(KP,3)*CCT
  1234.       CYZ(KP)=UMI(KP,2)*UMI(KP,3)*CCT
  1235.  
  1236.  8018 CONTINUE
  1237.  
  1238.       ELSEIF(IDCENN.EQ.4)THEN
  1239.       DT19=DTM1*0.5D0
  1240.       DO 8009 K=KDEB,KFIN
  1241.       KP=K-KDEB+1
  1242.       XMB(KP)=XMH(KP)
  1243.       CXT(KP)=UMI(KP,1)*UMI(KP,1)*DT19
  1244.       CYT(KP)=UMI(KP,2)*UMI(KP,2)*DT19
  1245.       CZT(KP)=UMI(KP,3)*UMI(KP,3)*DT19
  1246.       CXY(KP)=UMI(KP,1)*UMI(KP,2)*DT19
  1247.       CXZ(KP)=UMI(KP,1)*UMI(KP,3)*DT19
  1248.       CYZ(KP)=UMI(KP,2)*UMI(KP,3)*DT19
  1249.  8009 CONTINUE
  1250.  
  1251.       ENDIF
  1252.  
  1253.       ELSEIF(IKOMP.EQ.1)THEN
  1254.  
  1255. C ALFA,AKSI,CCU,CCT utilisées seulement dans la boucle 8008
  1256.       IF(IDCENN.EQ.1)THEN
  1257.       DO 81081 K=KDEB,KFIN
  1258.       KP=K-KDEB+1
  1259.       XMB(KP)=XMH(KP)
  1260.       CXT(KP)=0.D0
  1261.       CYT(KP)=0.D0
  1262.       CZT(KP)=0.D0
  1263.       CXY(KP)=0.D0
  1264.       CXZ(KP)=0.D0
  1265.       CYZ(KP)=0.D0
  1266.  
  1267.       DXT(KP)=0.D0
  1268.       DYT(KP)=0.D0
  1269.       DZT(KP)=0.D0
  1270.       DXY(KP)=0.D0
  1271.       DXZ(KP)=0.D0
  1272.       DYZ(KP)=0.D0
  1273.  
  1274. 81081 CONTINUE
  1275.  
  1276.       ELSEIF(IDCENN.EQ.2)THEN
  1277.       DO 81008 K=KDEB,KFIN
  1278.       KP=K-KDEB+1
  1279.       HMK=2.D0*UM(KP)/BM(KP)
  1280.       XMB(KP)=HMK
  1281.       ALFA=UM(KP)*HMK/ALFT(KP)/2.D0
  1282.       AKSI=ALFA/3.D0
  1283.       IF(ALFA.GT.3.D0)AKSI=1.D0
  1284.       CCT=AKSI/BM(KP)
  1285.  
  1286.       HPK=2.D0*UP(KP)/BP(KP)
  1287.       ALFA=UP(KP)*HPK/ALFT(KP)/2.D0
  1288.       AKSI=ALFA/3.D0
  1289.       IF(ALFA.GT.3.D0)AKSI=1.D0
  1290.       CCP=AKSI/BP(KP)
  1291.       CPT=CCP-CCT
  1292.       CC2=0.D0
  1293.       IF(CPT.GE.0.D0)CC2=CPT
  1294.       CC2 = CC2*1.3D0
  1295.  
  1296.       CXT(KP)=UMI(KP,1)*CCT
  1297.       CYT(KP)=UMI(KP,2)*CCT
  1298.       CZT(KP)=UMI(KP,3)*CCT
  1299.       CXY(KP)=CCT
  1300.       CXZ(KP)=CCT
  1301.       CYZ(KP)=CCT
  1302.  
  1303.       DXT(KP)=UPI(KP,1)*UPI(KP,1)*CC2
  1304.       DYT(KP)=UPI(KP,2)*UPI(KP,2)*CC2
  1305.       DZT(KP)=UPI(KP,3)*UPI(KP,3)*CC2
  1306.       DXY(KP)=UPI(KP,1)*UPI(KP,2)*CC2
  1307.       DXZ(KP)=UPI(KP,1)*UPI(KP,3)*CC2
  1308.       DYZ(KP)=UPI(KP,2)*UPI(KP,3)*CC2
  1309.  
  1310. 81008 CONTINUE
  1311.  
  1312.       ELSEIF(IDCENN.EQ.3)THEN
  1313.       DO 81018 K=KDEB,KFIN
  1314.       KP=K-KDEB+1
  1315.       HMK=2.D0*UM(KP)/BM(KP)
  1316.       XMB(KP)=HMK
  1317.       ALFA=UM(KP)*HMK/ALFT(KP)/2.D0
  1318.       AKSI=ALFA/3.D0
  1319.       IF(ALFA.GT.3.D0)AKSI=1.D0
  1320.       CCT=AKSI/BM(KP)
  1321.  
  1322.       CXT(KP)=UMI(KP,1)*CCT
  1323.       CYT(KP)=UMI(KP,2)*CCT
  1324.       CZT(KP)=UMI(KP,3)*CCT
  1325.       CXY(KP)=CCT
  1326.       CXZ(KP)=CCT
  1327.       CYZ(KP)=CCT
  1328.  
  1329.       DXT(KP)=0.D0
  1330.       DYT(KP)=0.D0
  1331.       DZT(KP)=0.D0
  1332.       DXY(KP)=0.D0
  1333.       DXZ(KP)=0.D0
  1334.       DYZ(KP)=0.D0
  1335.  
  1336. 81018 CONTINUE
  1337.  
  1338.       ELSEIF(IDCENN.EQ.4)THEN
  1339.       DT19=DTM1*0.5D0
  1340.       DO 81009 K=KDEB,KFIN
  1341.       KP=K-KDEB+1
  1342.       XMB(KP)=XMH(KP)
  1343.  
  1344.       CXT(KP)=UMI(KP,1)*DT19
  1345.       CYT(KP)=UMI(KP,2)*DT19
  1346.       CZT(KP)=UMI(KP,3)*DT19
  1347.       CXY(KP)=DT19
  1348.       CXZ(KP)=DT19
  1349.       CYZ(KP)=DT19
  1350.  
  1351.       DXT(KP)=0.D0
  1352.       DYT(KP)=0.D0
  1353.       DZT(KP)=0.D0
  1354.       DXY(KP)=0.D0
  1355.       DXZ(KP)=0.D0
  1356.       DYZ(KP)=0.D0
  1357.  
  1358. 81009 CONTINUE
  1359.  
  1360.       ENDIF
  1361.  
  1362.       ENDIF
  1363.  
  1364. C           ***********************
  1365. C
  1366. C     AVANT CALCUL DT
  1367. C
  1368.  
  1369.       IF(IKOMP.EQ.0)THEN
  1370.  
  1371. C*IBMDIR* PREFER SCALAR
  1372.       DO 8010 K=KDEB,KFIN
  1373.       KP=K-KDEB+1
  1374.       DT0=DT
  1375.  
  1376.       DT1=2.D0/
  1377.      & ( (UMI(KP,1)*UMI(KP,1))/(ALFT(KP)+CXT(KP))
  1378.      & + (UMI(KP,2)*UMI(KP,2))/(ALFT(KP)+CYT(KP))
  1379.      & + (UMI(KP,3)*UMI(KP,3))/(ALFT(KP)+CZT(KP)) )
  1380.  
  1381.       DT2=0.5/
  1382.      & ( (ALFT(KP)+CXT(KP))*AL2(KP)
  1383.      &  +(ALFT(KP)+CYT(KP))*AH2(KP)
  1384.      &  +(ALFT(KP)+CZT(KP))*AP2(KP) )
  1385.  
  1386.       IF(DT1.LT.DT)DT=DT1
  1387.       IF(DT2.LT.DT)DT=DT2
  1388. C     IF(DT3.LT.DT)DT=DT3
  1389.       IF(DT.NE.DT0) THEN
  1390.          DTT1=DT1
  1391.          DTT2=DT2
  1392.          DTT3=0.D0
  1393.          DIAEL=XMB(KP)
  1394.          NUEL=K
  1395.       END IF
  1396.  8010 CONTINUE
  1397.  
  1398.       ELSEIF(IKOMP.EQ.1)THEN
  1399.  
  1400. C*IBMDIR* PREFER SCALAR
  1401.       DO 8110 K=KDEB,KFIN
  1402.       KP=K-KDEB+1
  1403.       DT0=DT
  1404.  
  1405.       DT1=2.D0/
  1406.      & ( (UMI(KP,1)*UMI(KP,1))/(ALFT(KP)+CXT(KP)*UMI(KP,1)+DXT(KP))
  1407.      & + (UMI(KP,2)*UMI(KP,2))/(ALFT(KP)+CYT(KP)*UMI(KP,2)+DYT(KP))
  1408.      & + (UMI(KP,3)*UMI(KP,3))/(ALFT(KP)+CZT(KP)*UMI(KP,3)+DZT(KP)) )
  1409.  
  1410.       DT2=0.5/
  1411.      & ( (ALFT(KP)+CXT(KP)*UMI(KP,1)+DXT(KP))*AL2(KP)
  1412.      &  +(ALFT(KP)+CYT(KP)*UMI(KP,2)+DYT(KP))*AH2(KP)
  1413.      &  +(ALFT(KP)+CZT(KP)*UMI(KP,3)+DZT(KP))*AP2(KP) )
  1414.  
  1415.       IF(DT1.LT.DT)DT=DT1
  1416.       IF(DT2.LT.DT)DT=DT2
  1417. C     IF(DT3.LT.DT)DT=DT3
  1418.       IF(DT.NE.DT0) THEN
  1419.          DTT1=DT1
  1420.          DTT2=DT2
  1421.          DTT3=0.D0
  1422.          DIAEL=XMB(KP)
  1423.          NUEL=K
  1424.       ENDIF
  1425.  8110 CONTINUE
  1426.  
  1427.       ENDIF
  1428.  
  1429. C Le coeur du calcul ...
  1430.  
  1431.       IF(IKOMP.EQ.0)THEN
  1432.  
  1433.       DO 81094 I=1,NP
  1434.       DO 81093 J= 1,NP
  1435. C*IBMDIR* PREFER VECTOR
  1436.       DO 8014 K=KDEB,KFIN
  1437.       KP=K-KDEB+1
  1438. C     GEO1=(CF1(KP)*Q1(J,I)+CF2(KP)*Q2(J,I)+CF3(KP)*Q3(J,I))*CUB8
  1439.       GEO1=CFM(KP)*QGGT(J,I)*CUB8
  1440.       ZVGG=AIRE(KP)*(
  1441.      &     HR(K,I,1)*HR(K,J,1)*CXT(KP)
  1442.      &+    HR(K,I,2)*HR(K,J,2)*CYT(KP)
  1443.      &+    HR(K,I,3)*HR(K,J,3)*CZT(KP)
  1444.      &+    HR(K,I,1)*HR(K,J,2)*CXY(KP)
  1445.      &+    HR(K,I,2)*HR(K,J,1)*CXY(KP)
  1446.      &+    HR(K,I,1)*HR(K,J,3)*CXZ(KP)
  1447.      &+    HR(K,I,3)*HR(K,J,1)*CXZ(KP)
  1448.      &+    HR(K,I,2)*HR(K,J,3)*CYZ(KP)
  1449.      &+    HR(K,I,3)*HR(K,J,2)*CYZ(KP) )
  1450.      &+    (CXT(KP)+CYT(KP)+CZT(KP))/3.D0*GEO1
  1451. C    &+    (CXT(KP)+CYT(KP)+CZT(KP))/3.D0*XMH(KP)*QGGT(J,I)*CUB8
  1452.  
  1453.       ZVGT=AIRE(KP)*(
  1454.      &     HR(K,I,1)*HR(K,J,1)*ALFT(KP)
  1455.      &+    HR(K,I,2)*HR(K,J,2)*ALFT(KP)
  1456.      &+    HR(K,I,3)*HR(K,J,3)*ALFT(KP) )
  1457.      &+    ALFT(KP)*XMH(KP)*GEO1
  1458. C    &+    ALFT(KP)*XMH(KP)*QGGT(J,I)*CUB8
  1459.  
  1460.       V2=UMI(KP,1)*DR(KP,I)*HR(K,J,1)+UMI(KP,2)*DR(KP,I)*HR(K,J,2)
  1461.      &  +UMI(KP,3)*DR(KP,I)*HR(K,J,3)
  1462.  
  1463.       SBF(KP,I)=SBF(KP,I)+TETAC(KP,J)*(ZVGG+V2)+ TETAD(KP,J)*ZVGT
  1464.  
  1465.  8014 CONTINUE
  1466. 81093 CONTINUE
  1467. 81094 CONTINUE
  1468.  
  1469.       ELSEIF(IKOMP.EQ.1)THEN
  1470.  
  1471.       DO 81096 I=1,NP
  1472.       DO 81095 J= 1,NP
  1473. C*IBMDIR* PREFER VECTOR
  1474.       DO 8015 K=KDEB,KFIN
  1475.       KP=K-KDEB+1
  1476. C     GEO1=(CF1(KP)*Q1(J,I)+CF2(KP)*Q2(J,I)+CF3(KP)*Q3(J,I))*CUB8
  1477.       GEO1=CFM(KP)*QGGT(J,I)*CUB8
  1478.       ZVGG=AIRE(KP)*(
  1479.      &     HR(K,I,1)*HR(K,J,1)*(CXT(KP)*UIX(KP,J)+DXT(KP))
  1480.      &+    HR(K,I,2)*HR(K,J,2)*(CYT(KP)*UIY(KP,J)+DYT(KP))
  1481.      &+    HR(K,I,3)*HR(K,J,3)*(CZT(KP)*UIZ(KP,J)+DZT(KP))
  1482.      &+    HR(K,I,1)*HR(K,J,2)*(CXY(KP)*UMI(KP,1)*UIY(KP,J)+DXY(KP))
  1483.      &+    HR(K,I,2)*HR(K,J,1)*(CXY(KP)*UMI(KP,2)*UIX(KP,J)+DXY(KP))
  1484.      &+    HR(K,I,1)*HR(K,J,3)*(CXZ(KP)*UMI(KP,1)*UIZ(KP,J)+DXZ(KP))
  1485.      &+    HR(K,I,3)*HR(K,J,1)*(CXZ(KP)*UMI(KP,3)*UIX(KP,J)+DXZ(KP))
  1486.      &+    HR(K,I,2)*HR(K,J,3)*(CYZ(KP)*UMI(KP,2)*UIZ(KP,J)+DYZ(KP))
  1487.      &+    HR(K,I,3)*HR(K,J,2)*(CYZ(KP)*UMI(KP,3)*UIY(KP,J)+DYZ(KP))
  1488.      &+    (CXT(KP)*UIX(KP,J)+DXT(KP)+CYT(KP)*UIY(KP,J)+DYT(KP)
  1489.      &+     CZT(KP)*UIZ(KP,J)))/3.D0*GEO1
  1490. C    &+    (CXT(KP)+CYT(KP)+CZT(KP))/3.D0*XMH(KP)*QGGT(J,I)*CUB8
  1491.  
  1492.       ZVGT=AIRE(KP)*(
  1493.      &     HR(K,I,1)*HR(K,J,1)*ALFT(KP)
  1494.      &+    HR(K,I,2)*HR(K,J,2)*ALFT(KP)
  1495.      &+    HR(K,I,3)*HR(K,J,3)*ALFT(KP) )
  1496.      &+    ALFT(KP)*XMH(KP)*GEO1
  1497. C    &+    ALFT(KP)*XMH(KP)*QGGT(J,I)*CUB8
  1498.  
  1499.       V2=(UIX(KP,J)*HR(K,J,1)+UIY(KP,J)*HR(K,J,2)+
  1500.      & UIZ(KP,J)*HR(K,J,3))*DR(KP,I)
  1501.  
  1502.       SBF(KP,I)=SBF(KP,I)+TETAC(KP,J)*(ZVGG+V2)+ TETAD(KP,J)*ZVGT
  1503.  
  1504.  8015 CONTINUE
  1505. 81095 CONTINUE
  1506. 81096 CONTINUE
  1507.  
  1508.       ENDIF
  1509.  
  1510. C
  1511. C     Fin de l'accumulation dans SBF.
  1512. C     On ajoute ces incréments G.
  1513. C
  1514.       DO 81097 I=1,NP
  1515. C*IBMDIR* PREFER VECTOR
  1516.       DO 8017 K=KDEB,KFIN
  1517.       KP=K-KDEB+1
  1518.       NF=IPADS(LE(I,K))
  1519.       G(NF) = G(NF)-SBF(KP,I)
  1520.  8017 CONTINUE
  1521. 81097 CONTINUE
  1522.  
  1523.  8001 CONTINUE
  1524.  
  1525. C     WRITE(6,*)' G DANS YCTSCL '
  1526. C     WRITE(6,1984)(M,G(M),M=1,NPTD)
  1527.  
  1528. C     CALL ARRET(0)
  1529.       IPAS=1
  1530.       RETURN
  1531.  1002 FORMAT(10(1X,1PE11.4))
  1532.       END
  1533.  
  1534.  
  1535.  
  1536.  
  1537.  
  1538.  
  1539.  
  1540.  
  1541.  
  1542.  

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