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

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