Télécharger la plaquette Cast3M


fr en

Télécharger zjohns.eso

Retour à la liste

Numérotation des lignes : 

zjohns

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

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