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

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