Télécharger la plaquette Cast3M


fr en

Télécharger yctsca.eso

Retour à la liste

Numérotation des lignes : 

yctsca

  1. C YCTSCA    SOURCE    CHAT      05/01/13    04:16:28     5004
  2.       SUBROUTINE YCTSCA
  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. C! Attention : IAXI et RPG ne sont pas utilisées...
  12.      &     (HR,RPG,DRR,LE,NEL,K0,NPTD,IES,NP,IAXI,
  13.      &     IPADL,IKOMP,IKAS,IPADS,IPADD,
  14.      &     ALFE,IND1,UN,INDU,NPTS,TN,QE,IKS,
  15.      &     HRN,G,ALT,SGT,
  16.      &     VOLU,COTE,NELZ,IDCEN,
  17.      &     DTM1,DT,DTT1,DTT2,NUEL,DIAEL)
  18.  
  19.       IMPLICIT INTEGER(I-N)
  20.       IMPLICIT REAL*8 (A-H,O-Z)
  21.  
  22. C***********************************************************************
  23. C
  24. C  CE SP DISCRETISE UNE EQUATION GENERALE DE TRANSPORT-DIFFUSION AVEC
  25. C  SOURCE.
  26. C    EN 2D SUR LES ELEMENTS QUA4 ET TRI3      PLAN OU AXI
  27. C    EN 3D SUR LES ELEMENTS CUB8 ET PRI6
  28. C  LES OPERATEURS SONT "SOUS-INTEGRES"
  29. C
  30. C
  31. C                        APPELE PAR YTSCAL
  32. C
  33. C
  34. C***********************************************************************
  35.  
  36. C
  37. -INC PPARAM
  38. -INC CCOPTIO
  39. -INC CCVQUA4
  40. -INC CCREEL
  41. C
  42. C Longueur des registres vectoriels de la machine cible
  43. C On prend 64 pour ne pas augmenter la taille des tableaux
  44. C nécessaires à la vectorisation.
  45. C
  46.       PARAMETER(LRV=64)
  47.  
  48.       DIMENSION UN(NPTS,IES),HRN(NPTD),TN(NPTD)
  49.       DIMENSION COTE(NELZ,IES),VOLU(NELZ),QE(*)
  50.       DIMENSION ALFE(*),ALT(*)
  51.  
  52.       DIMENSION IPADL(*),LE(NP,*),IPADS(*),IPADD(*)
  53.       DIMENSION HR(NEL,NP,IES),RPG(*),DRR(NP,NEL)
  54.  
  55.       DIMENSION QGGT(8,8),Q1(8,8),Q2(8,8),Q3(8,8)
  56.  
  57.       DIMENSION AIRE(LRV)
  58.       DIMENSION AL2(LRV),AH2(LRV),AP2(LRV)
  59.       DIMENSION AL(LRV),AH(LRV),AP(LRV)
  60.       DIMENSION ALFT(LRV),QT(LRV)
  61.       DIMENSION XMB(LRV),XMH(LRV)
  62.       DIMENSION CFM(LRV)
  63. C     DIMENSION CF1(LRV),CF2(LRV),CF3(LRV)
  64.       DIMENSION DR(LRV,9)
  65.       DIMENSION UM(LRV),UP(LRV)
  66.       DIMENSION UIX(LRV,9),UIY(LRV,9),UIZ(LRV,9)
  67.       DIMENSION TETAC(LRV,9),TETAD(LRV,9)
  68.       DIMENSION UMI(LRV,3),UPI(LRV,3)
  69.       DIMENSION CXT(LRV),CYT(LRV),CXY(LRV)
  70.       DIMENSION DXT(LRV),DYT(LRV),DXY(LRV)
  71.       DIMENSION CZT(LRV),CXZ(LRV),CYZ(LRV)
  72.       DIMENSION DZT(LRV),DXZ(LRV),DYZ(LRV)
  73.       DIMENSION SBF(LRV,9)
  74.       DIMENSION GRADT(LRV,3)
  75.       DIMENSION BM(LRV),BP(LRV)
  76.  
  77.       REAL*8    G(*)
  78.  
  79.       SAVE IPAS,QGGT,Q1,Q2,Q3
  80.  
  81.       DATA CD/1.D0/
  82.  
  83.       DATA IPAS/0/
  84. C************************************************************************
  85. C
  86. C           INITIALISATIONS DIVERSES
  87. C
  88. C     WRITE(IOIMP,*)' debut YCTSCL',IES
  89.       ZERMA=1.D-20
  90.  
  91.       NK=K0
  92. C                           ********
  93. C                           *  2D  *
  94. C                           ********
  95.  
  96.       IF(IES.EQ.3)GO TO 10
  97.  
  98. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  99.  
  100. C                           DIFFERENCES TRIANGLE / QUADRANGLE
  101.       QUA4=0.D0
  102.       IF(NP.EQ.4)QUA4=1.D0
  103. C
  104. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  105.  
  106. C
  107. C Calcul du nombre de paquets de LRV éléments
  108. C
  109.       NNN=MOD(NEL,LRV)
  110.       IF(NNN.EQ.0) NPACK=NEL/LRV
  111.       IF(NNN.NE.0) NPACK=1+(NEL-NNN)/LRV
  112.       KPACKD=1
  113.       KPACKF=NPACK
  114. C
  115. C ******* BOUCLE SUR LES PAQUETS DE LRV ELEMENTS **********
  116. C
  117.       DO 7001 KPACK=KPACKD,KPACKF
  118. C
  119. C ======= A L'INTERIEUR DE CHAQUE PAQUET DE LRV ELEMENTS =======
  120. C
  121. C 1.D0 Calcul des limites du paquet courant.
  122.          KDEB=1+(KPACK-1)*LRV
  123.          KFIN=MIN(NEL,KDEB+LRV-1)
  124. C
  125.          DO 7002 K=KDEB,KFIN
  126.             KP=K-KDEB+1
  127.             NK=K+K0
  128.             NK1=(1-IND1)*(NK-1)+1
  129.             ALFT(KP)=ALFE(NK1)+ZERMA
  130.             AIRE(KP)=VOLU(NK)
  131.             AL(KP)=COTE(NK,1)
  132.             AH(KP)=COTE(NK,2)
  133.             AL2(KP)=1.D0/AL(KP)/AL(KP)
  134.             AH2(KP)=1.D0/AH(KP)/AH(KP)
  135.             XMH(KP)=(AL(KP)+AH(KP))/2.D0
  136.  7002    CONTINUE
  137.  
  138.          IF((IKOMP.EQ.0.AND.IKAS.EQ.5).OR.
  139.      &        (IKOMP.EQ.1.AND.IKAS.EQ.6))THEN
  140.             DO 7003 K=KDEB,KFIN
  141.                KP=K-KDEB+1
  142.                NK=K+K0
  143.                ALFT(KP)=ALFT(KP)+ALT(NK)/SGT
  144.  7003       CONTINUE
  145.          ENDIF
  146.  
  147. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  148.  
  149. C Initialisation des UMI avant accumulation
  150.          DO 7005 K=KDEB,KFIN
  151.             KP=K-KDEB+1
  152.             UMI(KP,1)=ZERMA
  153.             UMI(KP,2)=ZERMA
  154.             UPI(KP,1)=ZERMA
  155.             UPI(KP,2)=ZERMA
  156.             GRADT(KP,1)=ZERMA
  157.             GRADT(KP,2)=ZERMA
  158.  7005    CONTINUE
  159.  
  160.          DO 7006 I=1,NP
  161. C*IBMDIR* PREFER VECTOR
  162.             DO 70061 K=KDEB,KFIN
  163.                KP=K-KDEB+1
  164.                NF=IPADL(LE(I,K))
  165.                NFD=IPADD(LE(I,K))
  166.                NFU=IPADS(LE(I,K))
  167.                NFU=(1-INDU)*(NFU-1)+1
  168.                DR(KP,I)=DRR(I,K)
  169.                UIX(KP,I)=UN(NFU,1)
  170.                UIY(KP,I)=UN(NFU,2)
  171.                UMI(KP,1)=UMI(KP,1)+UN(NFU,1)*DR(KP,I)
  172.                UMI(KP,2)=UMI(KP,2)+UN(NFU,2)*DR(KP,I)
  173.  
  174.                TETAC(KP,I)=HRN(NF)
  175.                TETAD(KP,I)=TN(NFD)
  176.                GRADT(KP,1)=GRADT(KP,1)+HR(K,I,1)*TETAC(KP,I)
  177.                GRADT(KP,2)=GRADT(KP,2)+HR(K,I,2)*TETAC(KP,I)
  178. 70061       CONTINUE
  179.  7006    CONTINUE
  180.  
  181. C
  182. C Initialisation de la variable d'accumulation SBF au terme source
  183. C
  184. C       WRITE(IOIMP,*)' IKomp,ikas=',IKomp,ikas
  185. C       WRITE(IOIMP,*)' IKS,IND1,INDU=',IKS,IND1,INDU
  186.  
  187.          IF(IKOMP.EQ.0)THEN
  188.  
  189.             DO 70021 K=KDEB,KFIN
  190.                KP=K-KDEB+1
  191.                NK=K+K0
  192.                NKS=(1-IKS)*(NK-1)+1
  193.                QT(KP)=QE(NKS)
  194. 70021       CONTINUE
  195.  
  196.             DO 70062 I=1,NP
  197. C*IBMDIR* PREFER VECTOR
  198.                DO 70162 K=KDEB,KFIN
  199.                   KP=K-KDEB+1
  200.                   SBF(KP,I)=-QT(KP)*DR(KP,I)
  201. 70162          CONTINUE
  202. 70062       CONTINUE
  203.  
  204.          ELSEIF(IKOMP.EQ.1)THEN
  205.  
  206.             DO 70023 K=KDEB,KFIN
  207.                KP=K-KDEB+1
  208.                NK=K+K0
  209.                NKS=(1-IKS)*(NK-1)+1
  210.                QT(KP)=QE(NKS)
  211. 70023       CONTINUE
  212.  
  213.             DO 70064 I=1,NP
  214. C*IBMDIR* PREFER VECTOR
  215.                DO 70164 K=KDEB,KFIN
  216.                   KP=K-KDEB+1
  217.                   SBF(KP,I)=-QT(KP)*DR(KP,I)
  218. 70164          CONTINUE
  219. 70064       CONTINUE
  220.  
  221.          ENDIF
  222.  
  223.          DO 7007 K=KDEB,KFIN
  224.             KP=K-KDEB+1
  225.             UMI(KP,1)=UMI(KP,1)/AIRE(KP)
  226.             UMI(KP,2)=UMI(KP,2)/AIRE(KP)
  227.             UM(KP)=UMI(KP,1)*UMI(KP,1)+UMI(KP,2)*UMI(KP,2)
  228.             UM(KP)=SQRT(UM(KP))
  229.             A=GRADT(KP,1)*GRADT(KP,1)+GRADT(KP,2)*GRADT(KP,2)
  230.             A=SQRT(A)+XPETIT
  231.             GRADT(KP,1)=GRADT(KP,1)/A
  232.             GRADT(KP,2)=GRADT(KP,2)/A
  233.             UP(KP)=GRADT(KP,1)*UMI(KP,1)+GRADT(KP,2)*UMI(KP,2)
  234.             UPI(KP,1)=UP(KP)*GRADT(KP,1)
  235.             UPI(KP,2)=UP(KP)*GRADT(KP,2)
  236.  7007    CONTINUE
  237.  
  238.          DO 70074 K=KDEB,KFIN
  239.             KP=K-KDEB+1
  240.             BMX=UMI(KP,1)/AL(KP)
  241.             BMY=UMI(KP,2)/AH(KP)
  242.             BM(KP)=BMX*BMX+BMY*BMY
  243.             BM(KP)=SQRT(BM(KP))+XPETIT
  244.             BPX=UPI(KP,1)/AL(KP)
  245.             BPY=UPI(KP,2)/AH(KP)
  246.             BP(KP)=BPX*BPX+BPY*BPY
  247.             BP(KP)=SQRT(BP(KP))+XPETIT
  248. 70074    CONTINUE
  249.  
  250. C
  251. C     AVANT DECENTREMENT
  252. C
  253.  
  254.          IF(IKOMP.EQ.0)THEN
  255.  
  256. C ALFA,AKSI,CCU,CCT utilisées seulement dans la boucle 7008
  257.             IF(IDCEN.EQ.1)THEN
  258.                DO 70081 K=KDEB,KFIN
  259.                   KP=K-KDEB+1
  260.                   XMB(KP)=XMH(KP)
  261.                   CXT(KP)=0.D0
  262.                   CYT(KP)=0.D0
  263.                   CXY(KP)=0.D0
  264. 70081          CONTINUE
  265.  
  266.             ELSEIF(IDCEN.EQ.2)THEN
  267.                DO 7008 K=KDEB,KFIN
  268.                   KP=K-KDEB+1
  269.                   HMK=2.D0*UM(KP)/BM(KP)
  270.                   XMB(KP)=HMK
  271.                   ALFA=UM(KP)*HMK/ALFT(KP)/2.D0
  272.                   AKSI=ALFA/3.D0
  273.                   IF(ALFA.GT.3.D0)AKSI=1.D0
  274.                   CCT=AKSI/BM(KP)
  275.  
  276.                   HPK=2.D0*UP(KP)/BP(KP)
  277.                   ALFA=UP(KP)*HPK/ALFT(KP)/2.D0
  278.                   AKSI=ALFA/3.D0
  279.                   IF(ALFA.GT.3.D0)AKSI=1.D0
  280.                   CCP=AKSI/BP(KP)
  281.                   CPT=CCP-CCT
  282.                   CC2=0.D0
  283.                   IF(CPT.GE.0.D0)CC2=CPT
  284.  
  285.                   CXT(KP)=(UMI(KP,1)*UMI(KP,1)*CCT
  286.      $                    +UPI(KP,1)*UPI(KP,1)*CC2)
  287.                   CYT(KP)=(UMI(KP,2)*UMI(KP,2)*CCT
  288.      $                    +UPI(KP,2)*UPI(KP,2)*CC2)
  289.                   CXY(KP)=(UMI(KP,1)*UMI(KP,2)*CCT
  290.      $                    +UPI(KP,1)*UPI(KP,2)*CC2)
  291.  
  292.  7008          CONTINUE
  293.  
  294.             ELSEIF(IDCEN.EQ.3)THEN
  295.                DO 7018 K=KDEB,KFIN
  296.                   KP=K-KDEB+1
  297.                   HMK=2.D0*UM(KP)/BM(KP)
  298.                   XMB(KP)=HMK
  299.                   ALFA=UM(KP)*HMK/ALFT(KP)/2.D0
  300.                   AKSI=ALFA/3.D0
  301.                   IF(ALFA.GT.3.D0)AKSI=1.D0
  302.                   CCT=AKSI/BM(KP)
  303.  
  304.                   CXT(KP)=UMI(KP,1)*UMI(KP,1)*CCT
  305.                   CYT(KP)=UMI(KP,2)*UMI(KP,2)*CCT
  306.                   CXY(KP)=UMI(KP,1)*UMI(KP,2)*CCT
  307.  
  308.  7018          CONTINUE
  309.  
  310.             ELSEIF(IDCEN.EQ.4)THEN
  311.                DT19=DTM1*0.5D0
  312.                DO 7009 K=KDEB,KFIN
  313.                   KP=K-KDEB+1
  314.                   XMB(KP)=XMH(KP)
  315.                   CXT(KP)=UMI(KP,1)*UMI(KP,1)*DT19
  316.                   CYT(KP)=UMI(KP,2)*UMI(KP,2)*DT19
  317.                   CXY(KP)=UMI(KP,1)*UMI(KP,2)*DT19
  318.  7009          CONTINUE
  319.  
  320.             ENDIF
  321.  
  322.          ELSEIF(IKOMP.EQ.1)THEN
  323.  
  324. C ALFA,AKSI,CCU,CCT utilisées seulement dans la boucle 7008
  325.             IF(IDCEN.EQ.1)THEN
  326.                DO 71081 K=KDEB,KFIN
  327.                   KP=K-KDEB+1
  328.                   XMB(KP)=XMH(KP)
  329.  
  330.                   CXT(KP)=0.D0
  331.                   CYT(KP)=0.D0
  332.                   CXY(KP)=0.D0
  333.  
  334.                   DXT(KP)=0.D0
  335.                   DYT(KP)=0.D0
  336.                   DXY(KP)=0.D0
  337.  
  338. 71081          CONTINUE
  339.  
  340.             ELSEIF(IDCEN.EQ.2)THEN
  341.                DO 7108 K=KDEB,KFIN
  342.                   KP=K-KDEB+1
  343.                   HMK=2.D0*UM(KP)/BM(KP)
  344.                   XMB(KP)=HMK
  345.                   ALFA=UM(KP)*HMK/ALFT(KP)/2.D0
  346.                   AKSI=ALFA/3.D0
  347.                   IF(ALFA.GT.3.D0)AKSI=1.D0
  348.                   CCT=AKSI/BM(KP)
  349.  
  350.                   HPK=2.D0*UP(KP)/BP(KP)
  351.                   ALFA=UP(KP)*HPK/ALFT(KP)/2.D0
  352.                   AKSI=ALFA/3.D0
  353.                   IF(ALFA.GT.3.D0)AKSI=1.D0
  354.                   CCP=AKSI/BP(KP)
  355.                   CPT=CCP-CCT
  356.                   CC2=0.D0
  357.                   IF(CPT.GE.0.D0)CC2=CPT
  358.                   CC2 = CC2*1.3D0
  359.  
  360.                   CXT(KP)=UMI(KP,1)*CCT
  361.                   CYT(KP)=UMI(KP,2)*CCT
  362.                   CXY(KP)=CCT
  363.  
  364.                   DXT(KP)=UPI(KP,1)*UPI(KP,1)*CC2
  365.                   DYT(KP)=UPI(KP,2)*UPI(KP,2)*CC2
  366.                   DXY(KP)=UPI(KP,1)*UPI(KP,2)*CC2
  367.  
  368.  7108          CONTINUE
  369.  
  370.             ELSEIF(IDCEN.EQ.3)THEN
  371.                DO 71018 K=KDEB,KFIN
  372.                   KP=K-KDEB+1
  373.                   HMK=2.D0*UM(KP)/BM(KP)
  374.                   XMB(KP)=HMK
  375.                   ALFA=UM(KP)*HMK/ALFT(KP)/2.D0
  376.                   AKSI=ALFA/3.D0
  377.                   IF(ALFA.GT.3.D0)AKSI=1.D0
  378.                   CCT=AKSI/BM(KP)
  379.  
  380.                   CXT(KP)=UMI(KP,1)*CCT
  381.                   CYT(KP)=UMI(KP,2)*CCT
  382.                   CXY(KP)=CCT
  383.  
  384.                   DXT(KP)=0.D0
  385.                   DYT(KP)=0.D0
  386.                   DXY(KP)=0.D0
  387.  
  388. 71018          CONTINUE
  389.  
  390.             ELSEIF(IDCEN.EQ.4)THEN
  391.                DT19=DTM1*0.5D0
  392.                DO 71009 K=KDEB,KFIN
  393.                   KP=K-KDEB+1
  394.                   XMB(KP)=XMH(KP)
  395.                   CXT(KP)=UMI(KP,1)*DT19
  396.                   CYT(KP)=UMI(KP,2)*DT19
  397.                   CXY(KP)=DT19
  398.  
  399.                   DXT(KP)=0.D0
  400.                   DYT(KP)=0.D0
  401.                   DXY(KP)=0.D0
  402.  
  403. 71009          CONTINUE
  404.  
  405.             ENDIF
  406.  
  407.          ENDIF
  408. C           ***********************
  409. C
  410. C     AVANT CALCUL DT
  411. C
  412.          IF(IKOMP.EQ.0)THEN
  413.  
  414. C*IBMDIR* PREFER SCALAR
  415.             DO 7010 K=KDEB,KFIN
  416.                KP=K-KDEB+1
  417.                DT0=DT
  418.  
  419.                DT1=2.D0/
  420.      &              ( (UMI(KP,1)*UMI(KP,1))/(ALFT(KP)+CXT(KP))
  421.      &              + (UMI(KP,2)*UMI(KP,2))/(ALFT(KP)+CYT(KP)) )
  422.  
  423.                DT2=0.5D0/
  424.      &              ( (ALFT(KP)+CXT(KP))*AL2(KP)
  425.      &              +(ALFT(KP)+CYT(KP))*AH2(KP) )
  426.  
  427.                IF(DT1.LT.DT)DT=DT1
  428.                IF(DT2.LT.DT)DT=DT2
  429. C     IF(DT3.LT.DT)DT=DT3
  430.                IF(DT.NE.DT0) THEN
  431.                   DTT1=DT1
  432.                   DTT2=DT2
  433. C                  DTT3=0.D0
  434.                   DIAEL=XMB(KP)
  435.                   NUEL=K
  436.                END IF
  437.                CXY(KP)=CXY(KP)*CD
  438.  7010       CONTINUE
  439.  
  440.          ELSEIF(IKOMP.EQ.1)THEN
  441.  
  442. C*IBMDIR* PREFER SCALAR
  443.             DO 7110 K=KDEB,KFIN
  444.                KP=K-KDEB+1
  445.                DT0=DT
  446.  
  447.                DT1=2.D0/
  448.      &            ( (UMI(KP,1)*UMI(KP,1))
  449.      $             /(ALFT(KP)+CXT(KP)*UMI(KP,1)+DXT(KP))
  450.      $            + (UMI(KP,2)*UMI(KP,2))
  451.      $             /(ALFT(KP)+CYT(KP)*UMI(KP,2)+DYT(KP)) )
  452.  
  453.                DT2=0.5D0/
  454.      &              ( (ALFT(KP)+CXT(KP)*UMI(KP,1)+DXT(KP))*AL2(KP)
  455.      &              + (ALFT(KP)+CYT(KP)*UMI(KP,2)+DYT(KP))*AH2(KP) )
  456.  
  457.                IF(DT1.LT.DT)DT=DT1
  458.                IF(DT2.LT.DT)DT=DT2
  459. C     IF(DT3.LT.DT)DT=DT3
  460.                IF(DT.NE.DT0) THEN
  461.                   DTT1=DT1
  462.                   DTT2=DT2
  463. C                  DTT3=0.D0
  464.                   DIAEL=XMB(KP)
  465.                   NUEL=K
  466.                END IF
  467.                CXY(KP)=CXY(KP)*CD
  468.  7110       CONTINUE
  469.  
  470.          ENDIF
  471. C Le coeur du calcul ...
  472.  
  473.          IF(IKOMP.EQ.0)THEN
  474.  
  475.             DO 7014 I=1,NP
  476.                DO 70141 J= 1,NP
  477. C*IBMDIR* PREFER VECTOR
  478.                   DO 70142 K=KDEB,KFIN
  479.                      KP=K-KDEB+1
  480.                      ZVGG=AIRE(KP)*(
  481.      &                    HR(K,I,1)*HR(K,J,1)*CXT(KP)
  482.      &               +    HR(K,I,1)*HR(K,J,2)*CXY(KP)
  483.      &               +    HR(K,I,2)*HR(K,J,1)*CXY(KP)
  484.      &               +    HR(K,I,2)*HR(K,J,2)*CYT(KP)
  485.      &               +    (CXT(KP)*AL2(KP)+CYT(KP)*AH2(KP))/12.D0
  486.      &                    *VGGT(J,I)*QUA4  )
  487.  
  488. C    &+    HR(K,I,2)*HR(K,J,2)*CYT(KP) )
  489. C    &+    (CXT(KP)+CYT(KP))/24.D0*VGGT(J,I)*QUA4
  490.  
  491.                      ZVGT=AIRE(KP)*(
  492.      &                    HR(K,I,1)*HR(K,J,1)*ALFT(KP)
  493.      &               +    HR(K,I,2)*HR(K,J,2)*ALFT(KP)
  494.      &               +    (ALFT(KP)*AL2(KP)+ALFT(KP)*AH2(KP))/12.D0
  495.      &                    *VGGT(J,I)*QUA4  )
  496.  
  497. C    &+    HR(K,I,2)*HR(K,J,2)*ALFT(KP) )
  498. C    &+    ALFT(KP)/12.D0*VGGT(J,I)*QUA4
  499.  
  500. C?    V2=(UIX(KP,I)*HR(K,J,1)+UIY(KP,I)*HR(K,J,2))*DR(KP,I)
  501.                      V2=UMI(KP,1)*DR(KP,I)*HR(K,J,1)
  502.      $                 +UMI(KP,2)*DR(KP,I)*HR(K,J,2)
  503.  
  504.                      SBF(KP,I)=SBF(KP,I)+TETAC(KP,J)*(ZVGG+V2)
  505.      $                                  +TETAD(KP,J)*ZVGT
  506.  
  507. 70142             CONTINUE
  508. 70141          CONTINUE
  509.  7014       CONTINUE
  510.  
  511.          ELSEIF(IKOMP.EQ.1)THEN
  512.  
  513.             DO 7015 I=1,NP
  514.                DO 70151 J= 1,NP
  515. C*IBMDIR* PREFER VECTOR
  516.                   DO 70152 K=KDEB,KFIN
  517.                      KP=K-KDEB+1
  518.                      ZVGG=AIRE(KP)*(
  519.      &     HR(K,I,1)*HR(K,J,1)*(CXT(KP)*UIX(KP,J)+DXT(KP))
  520.      &+    HR(K,I,1)*HR(K,J,2)
  521.      &     *(CXY(KP)*UMI(KP,1)*UIY(KP,J)+DXY(KP))
  522.      &+    HR(K,I,2)*HR(K,J,1)
  523.      &     *(CXY(KP)*UMI(KP,2)*UIX(KP,J)+DXY(KP))
  524.      &+    HR(K,I,2)*HR(K,J,2)*(CYT(KP)*UIY(KP,J)+DYT(KP))
  525.      &+    ((CXT(KP)*UIX(KP,J)+DXT(KP))*AL2(KP)
  526.      &      +(CYT(KP)*UIY(KP,J)+DYT(KP))*AH2(KP))/12.D0
  527.      &     *VGGT(J,I)*QUA4  )
  528.  
  529. C    &+    HR(K,I,2)*HR(K,J,2)*CYT(KP) )
  530. C    &+    (CXT(KP)+CYT(KP))/24.D0*VGGT(J,I)*QUA4
  531.  
  532.                      ZVGT=AIRE(KP)*(
  533.      &     HR(K,I,1)*HR(K,J,1)*ALFT(KP)
  534.      &+    HR(K,I,2)*HR(K,J,2)*ALFT(KP)
  535.      &+    (ALFT(KP)*AL2(KP)+ALFT(KP)*AH2(KP))/12.D0
  536.      &     *VGGT(J,I)*QUA4  )
  537.  
  538. C    &+    HR(K,I,2)*HR(K,J,2)*ALFT(KP) )
  539. C    &+    ALFT(KP)/12.D0*VGGT(J,I)*QUA4
  540.  
  541.                      V2=(UIX(KP,J)*HR(K,J,1)
  542.      $                  +UIY(KP,J)*HR(K,J,2))*DR(KP,I)
  543.  
  544.                      SBF(KP,I)=SBF(KP,I)+TETAC(KP,J)*(ZVGG+V2)
  545.      $                                  +TETAD(KP,J)*ZVGT
  546.  
  547. 70152             CONTINUE
  548. 70151          CONTINUE
  549.  7015       CONTINUE
  550.  
  551.          ENDIF
  552. C
  553. C     Fin de l'accumulation dans SBF.
  554. C     On ajoute ces incréments G.
  555. C
  556.          DO 7017 I=1,NP
  557. C*IBMDIR* PREFER VECTOR
  558.             DO 70171 K=KDEB,KFIN
  559.                KP=K-KDEB+1
  560.                NF=IPADS(LE(I,K))
  561.                G(NF) = G(NF)-SBF(KP,I)
  562. 70171       CONTINUE
  563.  7017    CONTINUE
  564.  
  565.  7001 CONTINUE
  566.  
  567. C     WRITE(IOIMP,*)' G DANS YCTSCL '
  568. C     WRITE(IOIMP,1984)(M,G(M),M=1,NPTS)
  569.  1984 FORMAT(7(1X,I4,2X,1PE11.4))
  570.  
  571. C     CALL ARRET(0)
  572.       IPAS=1
  573.       RETURN
  574.  
  575. C                           ********
  576. C                           *  3D  *
  577. C                           ********
  578.  
  579.  10   CONTINUE
  580.  
  581.       IF(IPAS.EQ.0)CALL CALHRH(QGGT,Q1,Q2,Q3,IES)
  582.       CUB8=0.D0
  583.       IF(NP.EQ.8)CUB8=1.D0
  584.  
  585. C
  586. C Calcul du nombre de paquets de LRV éléments
  587. C
  588.       NNN=MOD(NEL,LRV)
  589.       IF(NNN.EQ.0) NPACK=NEL/LRV
  590.       IF(NNN.NE.0) NPACK=1+(NEL-NNN)/LRV
  591.       KPACKD=1
  592.       KPACKF=NPACK
  593. C
  594. C ******* BOUCLE SUR LES PAQUETS DE LRV ELEMENTS **********
  595. C
  596.       DO 8001 KPACK=KPACKD,KPACKF
  597. C
  598. C ======= A L'INTERIEUR DE CHAQUE PAQUET DE LRV ELEMENTS =======
  599. C
  600. C 1.D0 Calcul des limites du paquet courant.
  601.          KDEB=1+(KPACK-1)*LRV
  602.          KFIN=MIN(NEL,KDEB+LRV-1)
  603. C
  604.          DO 8002 K=KDEB,KFIN
  605.             KP=K-KDEB+1
  606.             NK=K+K0
  607.             NK1=(1-IND1)*(NK-1)+1
  608.             ALFT(KP)=ALFE(NK1)+ZERMA
  609.             AIRE(KP)=VOLU(NK)
  610.             AL(KP)=COTE(NK,1)
  611.             AH(KP)=COTE(NK,2)
  612.             AP(KP)=COTE(NK,3)
  613.             CFM(KP)=AL(KP)*AH(KP)/AP(KP)+AL(KP)*AP(KP)/AH(KP)+
  614.      &              AP(KP)*AH(KP)/AL(KP)
  615. C     CF1(KP)=AL(KP)*AH(KP)/AP(KP)
  616. C     CF2(KP)=AL(KP)*AP(KP)/AH(KP)
  617. C     CF3(KP)=AP(KP)*AH(KP)/AL(KP)
  618.             XMH(KP)=(AL(KP)+AH(KP)+AP(KP))/3.D0
  619.             AL2(KP)=1.D0/AL(KP)/AL(KP)
  620.             AH2(KP)=1.D0/AH(KP)/AH(KP)
  621.             AP2(KP)=1.D0/AP(KP)/AP(KP)
  622.  8002    CONTINUE
  623.          IF((IKOMP.EQ.0.AND.IKAS.EQ.5).OR.
  624.      &        (IKOMP.EQ.1.AND.IKAS.EQ.6))THEN
  625.             DO 8003 K=KDEB,KFIN
  626.                KP=K-KDEB+1
  627.                NK=K+K0
  628.                ALFT(KP)=ALFT(KP)+ALT(NK)/SGT
  629.  8003       CONTINUE
  630.          ENDIF
  631.  
  632. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  633.  
  634. C Initialisation des UMI avant accumulation
  635.          DO 8005 K=KDEB,KFIN
  636.             KP=K-KDEB+1
  637.             UMI(KP,1)=ZERMA
  638.             UMI(KP,2)=ZERMA
  639.             UMI(KP,3)=ZERMA
  640.             UPI(KP,1)=ZERMA
  641.             UPI(KP,2)=ZERMA
  642.             UPI(KP,3)=ZERMA
  643.             GRADT(KP,1)=ZERMA
  644.             GRADT(KP,2)=ZERMA
  645.             GRADT(KP,3)=ZERMA
  646.  8005    CONTINUE
  647.  
  648.          DO 8006 I=1,NP
  649. C*IBMDIR* PREFER VECTOR
  650.             DO 80061 K=KDEB,KFIN
  651.                KP=K-KDEB+1
  652.                NF=IPADL(LE(I,K))
  653.                NFD=IPADD(LE(I,K))
  654.                NFU=IPADS(LE(I,K))
  655.                NFU=(1-INDU)*(NFU-1)+1
  656.                DR(KP,I)=DRR(I,K)
  657.                UIX(KP,I)=UN(NFU,1)
  658.                UIY(KP,I)=UN(NFU,2)
  659.                UIZ(KP,I)=UN(NFU,3)
  660.                UMI(KP,1)=UMI(KP,1)+UN(NFU,1)*DR(KP,I)
  661.                UMI(KP,2)=UMI(KP,2)+UN(NFU,2)*DR(KP,I)
  662.                UMI(KP,3)=UMI(KP,3)+UN(NFU,3)*DR(KP,I)
  663.  
  664.                TETAC(KP,I)=HRN(NF)
  665.                TETAD(KP,I)=TN(NFD)
  666.                GRADT(KP,1)=GRADT(KP,1)+HR(K,I,1)*TETAC(KP,I)
  667.                GRADT(KP,2)=GRADT(KP,2)+HR(K,I,2)*TETAC(KP,I)
  668.                GRADT(KP,3)=GRADT(KP,3)+HR(K,I,3)*TETAC(KP,I)
  669. 80061       CONTINUE
  670.  8006    CONTINUE
  671.  
  672. C
  673. C Initialisation de la variable d'accumulation SBF au terme source
  674. C                                                                      M
  675. C       WRITE(IOIMP,*)' IKomp,ikas=',IKomp,ikas
  676. C       WRITE(IOIMP,*)' IKS,IND1,INDU=',IKS,IND1,INDU
  677.  
  678.          IF(IKOMP.EQ.0)THEN
  679.  
  680.             DO 80021 K=KDEB,KFIN
  681.                KP=K-KDEB+1
  682.                NK=K+K0
  683.                NKS=(1-IKS)*(NK-1)+1
  684.                QT(KP)=QE(NKS)
  685. 80021       CONTINUE
  686.  
  687.             DO 80062 I=1,NP
  688. C*IBMDIR* PREFER VECTOR
  689.                DO 80162 K=KDEB,KFIN
  690.                   KP=K-KDEB+1
  691.                   SBF(KP,I)=-QT(KP)*DR(KP,I)
  692. 80162          CONTINUE
  693. 80062       CONTINUE
  694.  
  695.          ELSEIF(IKOMP.EQ.1)THEN
  696.  
  697.             DO 80023 K=KDEB,KFIN
  698.                KP=K-KDEB+1
  699.                NK=K+K0
  700.                NKS=(1-IKS)*(NK-1)+1
  701.                QT(KP)=QE(NKS)
  702. 80023       CONTINUE
  703.  
  704.             DO 80064 I=1,NP
  705. C*IBMDIR* PREFER VECTOR
  706.                DO 80164 K=KDEB,KFIN
  707.                   KP=K-KDEB+1
  708.                   SBF(KP,I)=-QT(KP)*DR(KP,I)
  709. 80164          CONTINUE
  710. 80064       CONTINUE
  711.  
  712.          ENDIF
  713.  
  714.          DO 8007 K=KDEB,KFIN
  715.             KP=K-KDEB+1
  716.             UMI(KP,1)=UMI(KP,1)/AIRE(KP)
  717.             UMI(KP,2)=UMI(KP,2)/AIRE(KP)
  718.             UMI(KP,3)=UMI(KP,3)/AIRE(KP)
  719.             UM(KP)=UMI(KP,1)*UMI(KP,1)+UMI(KP,2)*UMI(KP,2)
  720.      &           +UMI(KP,3)*UMI(KP,3)
  721.             UM(KP)=SQRT(UM(KP))
  722.             A=GRADT(KP,1)*GRADT(KP,1)+GRADT(KP,2)*GRADT(KP,2)
  723.      &       +GRADT(KP,3)*GRADT(KP,3)
  724.             A=SQRT(A)+XPETIT
  725.             GRADT(KP,1)=GRADT(KP,1)/A
  726.             GRADT(KP,2)=GRADT(KP,2)/A
  727.             GRADT(KP,3)=GRADT(KP,3)/A
  728.             UP(KP)=GRADT(KP,1)*UMI(KP,1)+GRADT(KP,2)*UMI(KP,2)
  729.      &           +GRADT(KP,3)*UMI(KP,3)
  730.             UPI(KP,1)=UP(KP)*GRADT(KP,1)
  731.             UPI(KP,2)=UP(KP)*GRADT(KP,2)
  732.             UPI(KP,3)=UP(KP)*GRADT(KP,3)
  733.  8007    CONTINUE
  734.  
  735.          DO 80074 K=KDEB,KFIN
  736.             KP=K-KDEB+1
  737.             BMX=UMI(KP,1)/AL(KP)
  738.             BMY=UMI(KP,2)/AH(KP)
  739.             BMZ=UMI(KP,3)/AP(KP)
  740.             BM(KP)=BMX*BMX+BMY*BMY+BMZ*BMZ
  741.             BM(KP)=SQRT(BM(KP))+XPETIT
  742.             BPX=UPI(KP,1)/AL(KP)
  743.             BPY=UPI(KP,2)/AH(KP)
  744.             BPZ=UPI(KP,3)/AP(KP)
  745.             BP(KP)=BPX*BPX+BPY*BPY+BPZ*BPZ
  746.             BP(KP)=SQRT(BP(KP))+XPETIT
  747. 80074    CONTINUE
  748.  
  749. C
  750. C     AVANT DECENTREMENT
  751. C
  752.  
  753.          IF(IKOMP.EQ.0)THEN
  754.  
  755. C ALFA,AKSI,CCU,CCT utilisées seulement dans la boucle 8008
  756.             IF(IDCEN.EQ.1)THEN
  757.                DO 80081 K=KDEB,KFIN
  758.                   KP=K-KDEB+1
  759.                   XMB(KP)=XMH(KP)
  760.                   CXT(KP)=0.D0
  761.                   CYT(KP)=0.D0
  762.                   CZT(KP)=0.D0
  763.                   CXY(KP)=0.D0
  764.                   CXZ(KP)=0.D0
  765.                   CYZ(KP)=0.D0
  766. 80081          CONTINUE
  767.  
  768.             ELSEIF(IDCEN.EQ.2)THEN
  769.                DO 8008 K=KDEB,KFIN
  770.                   KP=K-KDEB+1
  771.                   HMK=2.D0*UM(KP)/BM(KP)
  772.                   XMB(KP)=HMK
  773.                   ALFA=UM(KP)*HMK/ALFT(KP)/2.D0
  774.                   AKSI=ALFA/3.D0
  775.                   IF(ALFA.GT.3.D0)AKSI=1.D0
  776.                   CCT=AKSI/BM(KP)
  777.  
  778.                   HPK=2.D0*UP(KP)/BP(KP)
  779.                   ALFA=UP(KP)*HPK/ALFT(KP)/2.D0
  780.                   AKSI=ALFA/3.D0
  781.                   IF(ALFA.GT.3.D0)AKSI=1.D0
  782.                   CCP=AKSI/BP(KP)
  783.                   CPT=CCP-CCT
  784.                   CC2=0.D0
  785.                   IF(CPT.GE.0.D0)CC2=CPT
  786.  
  787.                   CXT(KP)=(UMI(KP,1)*UMI(KP,1)*CCT
  788.      $                    +UPI(KP,1)*UPI(KP,1)*CC2)
  789.                   CYT(KP)=(UMI(KP,2)*UMI(KP,2)*CCT
  790.      $                    +UPI(KP,2)*UPI(KP,2)*CC2)
  791.                   CZT(KP)=(UMI(KP,3)*UMI(KP,3)*CCT
  792.      $                    +UPI(KP,3)*UPI(KP,3)*CC2)
  793.                   CXY(KP)=(UMI(KP,1)*UMI(KP,2)*CCT
  794.      $                    +UPI(KP,1)*UPI(KP,2)*CC2)
  795.                   CXZ(KP)=(UMI(KP,1)*UMI(KP,3)*CCT
  796.      $                    +UPI(KP,1)*UPI(KP,3)*CC2)
  797.                   CYZ(KP)=(UMI(KP,2)*UMI(KP,3)*CCT
  798.      $                    +UPI(KP,2)*UPI(KP,3)*CC2)
  799.  
  800.  8008          CONTINUE
  801.  
  802.             ELSEIF(IDCEN.EQ.3)THEN
  803.                DO 8018 K=KDEB,KFIN
  804.                   KP=K-KDEB+1
  805.                   HMK=2.D0*UM(KP)/BM(KP)
  806.                   XMB(KP)=HMK
  807.                   ALFA=UM(KP)*HMK/ALFT(KP)/2.D0
  808.                   AKSI=ALFA/3.D0
  809.                   IF(ALFA.GT.3.D0)AKSI=1.D0
  810.                   CCT=AKSI/BM(KP)
  811.  
  812.                   CXT(KP)=UMI(KP,1)*UMI(KP,1)*CCT
  813.                   CYT(KP)=UMI(KP,2)*UMI(KP,2)*CCT
  814.                   CZT(KP)=UMI(KP,3)*UMI(KP,3)*CCT
  815.                   CXY(KP)=UMI(KP,1)*UMI(KP,2)*CCT
  816.                   CXZ(KP)=UMI(KP,1)*UMI(KP,3)*CCT
  817.                   CYZ(KP)=UMI(KP,2)*UMI(KP,3)*CCT
  818.  
  819.  8018          CONTINUE
  820.  
  821.             ELSEIF(IDCEN.EQ.4)THEN
  822.                DT19=DTM1*0.5D0
  823.                DO 8009 K=KDEB,KFIN
  824.                   KP=K-KDEB+1
  825.                   XMB(KP)=XMH(KP)
  826.                   CXT(KP)=UMI(KP,1)*UMI(KP,1)*DT19
  827.                   CYT(KP)=UMI(KP,2)*UMI(KP,2)*DT19
  828.                   CZT(KP)=UMI(KP,3)*UMI(KP,3)*DT19
  829.                   CXY(KP)=UMI(KP,1)*UMI(KP,2)*DT19
  830.                   CXZ(KP)=UMI(KP,1)*UMI(KP,3)*DT19
  831.                   CYZ(KP)=UMI(KP,2)*UMI(KP,3)*DT19
  832.  8009          CONTINUE
  833.  
  834.             ENDIF
  835.  
  836.          ELSEIF(IKOMP.EQ.1)THEN
  837.  
  838. C ALFA,AKSI,CCU,CCT utilisées seulement dans la boucle 8008
  839.             IF(IDCEN.EQ.1)THEN
  840.                DO 81081 K=KDEB,KFIN
  841.                   KP=K-KDEB+1
  842.                   XMB(KP)=XMH(KP)
  843.                   CXT(KP)=0.D0
  844.                   CYT(KP)=0.D0
  845.                   CZT(KP)=0.D0
  846.                   CXY(KP)=0.D0
  847.                   CXZ(KP)=0.D0
  848.                   CYZ(KP)=0.D0
  849.  
  850.                   DXT(KP)=0.D0
  851.                   DYT(KP)=0.D0
  852.                   DZT(KP)=0.D0
  853.                   DXY(KP)=0.D0
  854.                   DXZ(KP)=0.D0
  855.                   DYZ(KP)=0.D0
  856.  
  857. 81081          CONTINUE
  858.  
  859.             ELSEIF(IDCEN.EQ.2)THEN
  860.                DO 81008 K=KDEB,KFIN
  861.                   KP=K-KDEB+1
  862.                   HMK=2.D0*UM(KP)/BM(KP)
  863.                   XMB(KP)=HMK
  864.                   ALFA=UM(KP)*HMK/ALFT(KP)/2.D0
  865.                   AKSI=ALFA/3.D0
  866.                   IF(ALFA.GT.3.D0)AKSI=1.D0
  867.                   CCT=AKSI/BM(KP)
  868.  
  869.                   HPK=2.D0*UP(KP)/BP(KP)
  870.                   ALFA=UP(KP)*HPK/ALFT(KP)/2.D0
  871.                   AKSI=ALFA/3.D0
  872.                   IF(ALFA.GT.3.D0)AKSI=1.D0
  873.                   CCP=AKSI/BP(KP)
  874.                   CPT=CCP-CCT
  875.                   CC2=0.D0
  876.                   IF(CPT.GE.0.D0)CC2=CPT
  877.                   CC2 = CC2*1.3D0
  878.  
  879.                   CXT(KP)=UMI(KP,1)*CCT
  880.                   CYT(KP)=UMI(KP,2)*CCT
  881.                   CZT(KP)=UMI(KP,3)*CCT
  882.                   CXY(KP)=CCT
  883.                   CXZ(KP)=CCT
  884.                   CYZ(KP)=CCT
  885.  
  886.                   DXT(KP)=UPI(KP,1)*UPI(KP,1)*CC2
  887.                   DYT(KP)=UPI(KP,2)*UPI(KP,2)*CC2
  888.                   DZT(KP)=UPI(KP,3)*UPI(KP,3)*CC2
  889.                   DXY(KP)=UPI(KP,1)*UPI(KP,2)*CC2
  890.                   DXZ(KP)=UPI(KP,1)*UPI(KP,3)*CC2
  891.                   DYZ(KP)=UPI(KP,2)*UPI(KP,3)*CC2
  892.  
  893. 81008          CONTINUE
  894.  
  895.             ELSEIF(IDCEN.EQ.3)THEN
  896.                DO 81018 K=KDEB,KFIN
  897.                   KP=K-KDEB+1
  898.                   HMK=2.D0*UM(KP)/BM(KP)
  899.                   XMB(KP)=HMK
  900.                   ALFA=UM(KP)*HMK/ALFT(KP)/2.D0
  901.                   AKSI=ALFA/3.D0
  902.                   IF(ALFA.GT.3.D0)AKSI=1.D0
  903.                   CCT=AKSI/BM(KP)
  904.  
  905.                   CXT(KP)=UMI(KP,1)*CCT
  906.                   CYT(KP)=UMI(KP,2)*CCT
  907.                   CZT(KP)=UMI(KP,3)*CCT
  908.                   CXY(KP)=CCT
  909.                   CXZ(KP)=CCT
  910.                   CYZ(KP)=CCT
  911.  
  912.                   DXT(KP)=0.D0
  913.                   DYT(KP)=0.D0
  914.                   DZT(KP)=0.D0
  915.                   DXY(KP)=0.D0
  916.                   DXZ(KP)=0.D0
  917.                   DYZ(KP)=0.D0
  918.  
  919. 81018          CONTINUE
  920.  
  921.             ELSEIF(IDCEN.EQ.4)THEN
  922.                DT19=DTM1*0.5D0
  923.                DO 81009 K=KDEB,KFIN
  924.                   KP=K-KDEB+1
  925.                   XMB(KP)=XMH(KP)
  926.  
  927.                   CXT(KP)=UMI(KP,1)*DT19
  928.                   CYT(KP)=UMI(KP,2)*DT19
  929.                   CZT(KP)=UMI(KP,3)*DT19
  930.                   CXY(KP)=DT19
  931.                   CXZ(KP)=DT19
  932.                   CYZ(KP)=DT19
  933.  
  934.                   DXT(KP)=0.D0
  935.                   DYT(KP)=0.D0
  936.                   DZT(KP)=0.D0
  937.                   DXY(KP)=0.D0
  938.                   DXZ(KP)=0.D0
  939.                   DYZ(KP)=0.D0
  940.  
  941. 81009          CONTINUE
  942.  
  943.             ENDIF
  944.  
  945.          ENDIF
  946.  
  947. C           ***********************
  948. C
  949. C     AVANT CALCUL DT
  950. C
  951.  
  952.          IF(IKOMP.EQ.0)THEN
  953.  
  954. C*IBMDIR* PREFER SCALAR
  955.             DO 8010 K=KDEB,KFIN
  956.                KP=K-KDEB+1
  957.                DT0=DT
  958.  
  959.                DT1=2.D0/
  960.      &              ( (UMI(KP,1)*UMI(KP,1))/(ALFT(KP)+CXT(KP))
  961.      &              + (UMI(KP,2)*UMI(KP,2))/(ALFT(KP)+CYT(KP))
  962.      &              + (UMI(KP,3)*UMI(KP,3))/(ALFT(KP)+CZT(KP)) )
  963.  
  964.                DT2=0.5D0/
  965.      &             ( (ALFT(KP)+CXT(KP))*AL2(KP)
  966.      &              +(ALFT(KP)+CYT(KP))*AH2(KP)
  967.      &              +(ALFT(KP)+CZT(KP))*AP2(KP) )
  968.  
  969.                IF(DT1.LT.DT)DT=DT1
  970.                IF(DT2.LT.DT)DT=DT2
  971. C     IF(DT3.LT.DT)DT=DT3
  972.                IF(DT.NE.DT0) THEN
  973.                   DTT1=DT1
  974.                   DTT2=DT2
  975. C                  DTT3=0.D0
  976.                   DIAEL=XMB(KP)
  977.                   NUEL=K
  978.                END IF
  979.  8010       CONTINUE
  980.  
  981.          ELSEIF(IKOMP.EQ.1)THEN
  982.  
  983. C*IBMDIR* PREFER SCALAR
  984.             DO 8110 K=KDEB,KFIN
  985.                KP=K-KDEB+1
  986.                DT0=DT
  987.  
  988.                DT1=2.D0/
  989.      &              ( (UMI(KP,1)*UMI(KP,1))
  990.      $               /(ALFT(KP)+CXT(KP)*UMI(KP,1)+DXT(KP))
  991.      $              + (UMI(KP,2)*UMI(KP,2))
  992.      $               /(ALFT(KP)+CYT(KP)*UMI(KP,2)+DYT(KP))
  993.      $              + (UMI(KP,3)*UMI(KP,3))
  994.      $               /(ALFT(KP)+CZT(KP)*UMI(KP,3)+DZT(KP)) )
  995.  
  996.                DT2=0.5D0/
  997.      &             ( (ALFT(KP)+CXT(KP)*UMI(KP,1)+DXT(KP))*AL2(KP)
  998.      &              +(ALFT(KP)+CYT(KP)*UMI(KP,2)+DYT(KP))*AH2(KP)
  999.      &              +(ALFT(KP)+CZT(KP)*UMI(KP,3)+DZT(KP))*AP2(KP) )
  1000.  
  1001.                IF(DT1.LT.DT)DT=DT1
  1002.                IF(DT2.LT.DT)DT=DT2
  1003. C     IF(DT3.LT.DT)DT=DT3
  1004.                IF(DT.NE.DT0) THEN
  1005.                   DTT1=DT1
  1006.                   DTT2=DT2
  1007. C                  DTT3=0.D0
  1008.                   DIAEL=XMB(KP)
  1009.                   NUEL=K
  1010.                ENDIF
  1011.  8110       CONTINUE
  1012.  
  1013.          ENDIF
  1014.  
  1015. C Le coeur du calcul ...
  1016.  
  1017.          IF(IKOMP.EQ.0)THEN
  1018.  
  1019.             DO 8014 I=1,NP
  1020.                DO 80141 J= 1,NP
  1021. C*IBMDIR* PREFER VECTOR
  1022.                   DO 80142 K=KDEB,KFIN
  1023.                      KP=K-KDEB+1
  1024. C     GEO1=(CF1(KP)*Q1(J,I)+CF2(KP)*Q2(J,I)+CF3(KP)*Q3(J,I))*CUB8
  1025.                      GEO1=CFM(KP)*QGGT(J,I)*CUB8
  1026.                      ZVGG=AIRE(KP)*(
  1027.      &                    HR(K,I,1)*HR(K,J,1)*CXT(KP)
  1028.      &                    +    HR(K,I,2)*HR(K,J,2)*CYT(KP)
  1029.      &                    +    HR(K,I,3)*HR(K,J,3)*CZT(KP)
  1030.      &                    +    HR(K,I,1)*HR(K,J,2)*CXY(KP)
  1031.      &                    +    HR(K,I,2)*HR(K,J,1)*CXY(KP)
  1032.      &                    +    HR(K,I,1)*HR(K,J,3)*CXZ(KP)
  1033.      &                    +    HR(K,I,3)*HR(K,J,1)*CXZ(KP)
  1034.      &                    +    HR(K,I,2)*HR(K,J,3)*CYZ(KP)
  1035.      &                    +    HR(K,I,3)*HR(K,J,2)*CYZ(KP) )
  1036.      &                    +    (CXT(KP)+CYT(KP)+CZT(KP))/3.D0*GEO1
  1037. C    &+    (CXT(KP)+CYT(KP)+CZT(KP))/3.D0*XMH(KP)*QGGT(J,I)*CUB8
  1038.  
  1039.                      ZVGT=AIRE(KP)*(
  1040.      &                    HR(K,I,1)*HR(K,J,1)*ALFT(KP)
  1041.      &                    +    HR(K,I,2)*HR(K,J,2)*ALFT(KP)
  1042.      &                    +    HR(K,I,3)*HR(K,J,3)*ALFT(KP) )
  1043.      &                    +    ALFT(KP)*XMH(KP)*GEO1
  1044. C    &+    ALFT(KP)*XMH(KP)*QGGT(J,I)*CUB8
  1045.  
  1046.                      V2=UMI(KP,1)*DR(KP,I)*HR(K,J,1)
  1047.      $                    +UMI(KP,2)*DR(KP,I)*HR(K,J,2)
  1048.      $                    +UMI(KP,3)*DR(KP,I)*HR(K,J,3)
  1049.  
  1050.                      SBF(KP,I)=SBF(KP,I)+TETAC(KP,J)*(ZVGG+V2)
  1051.      $                                  +TETAD(KP,J)*ZVGT
  1052.  
  1053. 80142             CONTINUE
  1054. 80141          CONTINUE
  1055.  8014       CONTINUE
  1056.  
  1057.          ELSEIF(IKOMP.EQ.1)THEN
  1058.  
  1059.             DO 8015 I=1,NP
  1060.                DO 80151 J= 1,NP
  1061. C*IBMDIR* PREFER VECTOR
  1062.                   DO 80152 K=KDEB,KFIN
  1063.                      KP=K-KDEB+1
  1064. C     GEO1=(CF1(KP)*Q1(J,I)+CF2(KP)*Q2(J,I)+CF3(KP)*Q3(J,I))*CUB8
  1065.                      GEO1=CFM(KP)*QGGT(J,I)*CUB8
  1066.                      ZVGG=AIRE(KP)*(
  1067.      &     HR(K,I,1)*HR(K,J,1)*(CXT(KP)*UIX(KP,J)+DXT(KP))
  1068.      &+    HR(K,I,2)*HR(K,J,2)*(CYT(KP)*UIY(KP,J)+DYT(KP))
  1069.      &+    HR(K,I,3)*HR(K,J,3)*(CZT(KP)*UIZ(KP,J)+DZT(KP))
  1070.      &+    HR(K,I,1)*HR(K,J,2)*(CXY(KP)*UMI(KP,1)*UIY(KP,J)+DXY(KP))
  1071.      &+    HR(K,I,2)*HR(K,J,1)*(CXY(KP)*UMI(KP,2)*UIX(KP,J)+DXY(KP))
  1072.      &+    HR(K,I,1)*HR(K,J,3)*(CXZ(KP)*UMI(KP,1)*UIZ(KP,J)+DXZ(KP))
  1073.      &+    HR(K,I,3)*HR(K,J,1)*(CXZ(KP)*UMI(KP,3)*UIX(KP,J)+DXZ(KP))
  1074.      &+    HR(K,I,2)*HR(K,J,3)*(CYZ(KP)*UMI(KP,2)*UIZ(KP,J)+DYZ(KP))
  1075.      &+    HR(K,I,3)*HR(K,J,2)*(CYZ(KP)*UMI(KP,3)*UIY(KP,J)+DYZ(KP))
  1076.      &+    (CXT(KP)*UIX(KP,J)+DXT(KP)+CYT(KP)*UIY(KP,J)+DYT(KP)
  1077.      &+     CZT(KP)*UIZ(KP,J)))/3.D0*GEO1
  1078. C    &+    (CXT(KP)+CYT(KP)+CZT(KP))/3.D0*XMH(KP)*QGGT(J,I)*CUB8
  1079.  
  1080.                      ZVGT=AIRE(KP)*(
  1081.      &     HR(K,I,1)*HR(K,J,1)*ALFT(KP)
  1082.      &+    HR(K,I,2)*HR(K,J,2)*ALFT(KP)
  1083.      &+    HR(K,I,3)*HR(K,J,3)*ALFT(KP) )
  1084.      &+    ALFT(KP)*XMH(KP)*GEO1
  1085. C    &+    ALFT(KP)*XMH(KP)*QGGT(J,I)*CUB8
  1086.  
  1087.                      V2=(UIX(KP,J)*HR(K,J,1)
  1088.      $                    +UIY(KP,J)*HR(K,J,2)
  1089.      $                    +UIZ(KP,J)*HR(K,J,3))*DR(KP,I)
  1090.  
  1091.                      SBF(KP,I)=SBF(KP,I)+TETAC(KP,J)*(ZVGG+V2)
  1092.      $                                  +TETAD(KP,J)*ZVGT
  1093.  
  1094. 80152             CONTINUE
  1095. 80151          CONTINUE
  1096.  8015       CONTINUE
  1097.  
  1098.          ENDIF
  1099.  
  1100. C
  1101. C     Fin de l'accumulation dans SBF.
  1102. C     On ajoute ces incréments G.
  1103. C
  1104.          DO 8017 I=1,NP
  1105. C*IBMDIR* PREFER VECTOR
  1106.             DO 80171 K=KDEB,KFIN
  1107.                KP=K-KDEB+1
  1108.                NF=IPADS(LE(I,K))
  1109.                G(NF) = G(NF)-SBF(KP,I)
  1110. 80171       CONTINUE
  1111.  8017    CONTINUE
  1112.  
  1113.  8001 CONTINUE
  1114.  
  1115. C     WRITE(IOIMP,*)' G DANS YCTSCL '
  1116. C     WRITE(IOIMP,1984)(M,G(M),M=1,NPTS)
  1117.  
  1118. C     CALL ARRET(0)
  1119.       IPAS=1
  1120.       RETURN
  1121.  1002 FORMAT(10(1X,1PE11.4))
  1122.       END
  1123.  
  1124.  
  1125.  
  1126.  
  1127.  
  1128.  
  1129.  

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