Télécharger la plaquette Cast3M


fr en

Télécharger ksupg1.eso

Retour à la liste

Numérotation des lignes : 

ksupg1

  1. C KSUPG1    SOURCE    CHAT      05/01/13    01:09:14     5004
  2.       SUBROUTINE KSUPG1(WTILDE,UMI,CHGLD,CHGLP,KDEB,KFIN,LRV,
  3.      &     TN,IPADI,UN,COEFT,NPTD,NEL,NP,DRR,HR,FN,
  4.      &     AIRE,AL,AH,AP,IDCEN,IPADU,LE,CHGL,IKOMP,
  5.      &     DTM1,DT,DTT1,DTT2,DIAEL,NUEL)
  6.       IMPLICIT INTEGER(I-N)
  7.       IMPLICIT REAL*8 (A-H,O-Z)
  8.  
  9. -INC PPARAM
  10. -INC CCOPTIO
  11. -INC CCREEL
  12.       DIMENSION IPADU(*),LE(NP,1),IPADI(*)
  13.       DIMENSION UN(NPTD,IDIM),TN(*)
  14.       DIMENSION AIRE(LRV),AL(LRV),AH(LRV),AP(LRV),WTILDE(LRV,NP)
  15.       DIMENSION COEFT(LRV),CHGLD(LRV),CHGLP(LRV)
  16.       DIMENSION HR(NEL,NP,IDIM),DRR(NP,NEL),FN(NP,1)
  17.       PARAMETER (LRV1=64)
  18.       DIMENSION UM(LRV1),UMI(LRV,3),UMJ(LRV1,3)
  19.       DIMENSION GRAD(LRV1,3)
  20.       DIMENSION UP(LRV1),UPI(LRV1,3)
  21.       DIMENSION BM(LRV1),BP(LRV1)
  22.       DIMENSION TO1(LRV1),TO2(LRV1)
  23.       DIMENSION CXT(LRV1),CYT(LRV1),CZT(LRV1)
  24.       DIMENSION CXY(LRV1),CXZ(LRV1),CYZ(LRV1),CFM(LRV1)
  25.       DIMENSION AL2(LRV1),AH2(LRV1),AP2(LRV1),XMB(LRV1)
  26.       DATA IKV/1/,NPG/1/
  27.       CC0=1.D0
  28.       xpeti2 =  sqrt(XPETIT)
  29.       IF(IKOMP.EQ.1)CC0=0.2D0
  30.       IF(LRV.NE.LRV1)STOP
  31.       DO 7005 N=1,IDIM
  32.          DO 70051 K=KDEB,KFIN
  33.             KP=K-KDEB+1
  34.             UMI(KP,N)=xpeti2
  35.             UMJ(KP,N)=xpeti2
  36.             UPI(KP,N)=xpeti2
  37.             GRAD(KP,N)=xpeti2
  38. 70051    CONTINUE
  39.  7005 CONTINUE
  40.       IF(IKV.EQ.1)THEN
  41.          DO 70060 N=1,IDIM
  42.             DO 70061 I=1,NP
  43. C*IBMDIR* PREFER VECTOR
  44.                DO 70062 K=KDEB,KFIN
  45.                   KP=K-KDEB+1
  46.                   NF=IPADU(LE(I,K))
  47.                   NG=IPADI(LE(I,K))
  48.                   UMI(KP,N)=UMI(KP,N)+UN(NF,N)*DRR(I,K)
  49.                   UMJ(KP,N)=UMJ(KP,N)+UN(NF,N)*DRR(I,K)
  50.                   GRAD(KP,N)=GRAD(KP,N)+HR(K,I,N)*TN(NG)
  51. 70062          CONTINUE
  52. 70061       CONTINUE
  53. 70060    CONTINUE
  54.          DO 7017 K=KDEB,KFIN
  55.             KP=K-KDEB+1
  56.             UMJ(KP,1)=UMJ(KP,1)/AIRE(KP)
  57.             UMJ(KP,2)=UMJ(KP,2)/AIRE(KP)
  58.             UM(KP)=UMJ(KP,1)*UMJ(KP,1)+UMJ(KP,2)*UMJ(KP,2)
  59.             UM(KP)=SQRT(UM(KP))
  60.  7017    CONTINUE
  61.       ELSE
  62.          DO 70160 N=1,IDIM
  63.             DO 70161 I=1,NP
  64. C*IBMDIR* PREFER VECTOR
  65.                DO 70162 K=KDEB,KFIN
  66.                   KP=K-KDEB+1
  67.                   NF=IPADU(LE(I,K))
  68.                   NG=IPADI(LE(I,K))
  69.                   UMI(KP,N)=UMI(KP,N)+UN(NF,N)*DRR(I,K)
  70.                   GRAD(KP,N)=GRAD(KP,N)+HR(K,I,N)*TN(NG)
  71.                   U=0.D0
  72.                   DO 7026 L=1,NPG
  73.                      U=U+UN(NF,N)*FN(I,L)
  74.  7026             CONTINUE
  75.                   UMJ(KP,N)=UMJ(KP,N)+U
  76. 70162          CONTINUE
  77. 70161       CONTINUE
  78. 70160    CONTINUE
  79.          DO 7027 K=KDEB,KFIN
  80.             KP=K-KDEB+1
  81.             UM(KP)=UMJ(KP,1)*UMJ(KP,1)+UMJ(KP,2)*UMJ(KP,2)
  82.             UM(KP)=SQRT(UM(KP))
  83.  7027    CONTINUE
  84.       ENDIF
  85. C                 ***********
  86. C                 *   2 D   *
  87. C                 ***********
  88.       IF(IDIM.EQ.2)THEN
  89.          QUA4=CHGL
  90.          DO 7004 K=KDEB,KFIN
  91.             KP=K-KDEB+1
  92.             AL2(KP)=1.D0/AL(KP)/AL(KP)
  93.             AH2(KP)=1.D0/AH(KP)/AH(KP)
  94.             XMB(KP)=(AL(KP)+AH(KP))/2.D0
  95.  7004    CONTINUE
  96.          DO 7007 K=KDEB,KFIN
  97.             KP=K-KDEB+1
  98.             UMI(KP,1)=UMI(KP,1)/AIRE(KP)
  99.             UMI(KP,2)=UMI(KP,2)/AIRE(KP)
  100.             AX=GRAD(KP,1)*GRAD(KP,1)+GRAD(KP,2)*GRAD(KP,2)
  101.             AX=SQRT(AX)+xpeti2
  102.             GRAD(KP,1)=GRAD(KP,1)/AX
  103.             GRAD(KP,2)=GRAD(KP,2)/AX
  104.             UP(KP)=(GRAD(KP,1)*UMJ(KP,1)+GRAD(KP,2)*UMJ(KP,2))
  105.             UPI(KP,1)=UP(KP)*GRAD(KP,1)
  106.             UPI(KP,2)=UP(KP)*GRAD(KP,2)
  107.  7007    CONTINUE
  108.          DO 70074 K=KDEB,KFIN
  109.             KP=K-KDEB+1
  110.             BMX=UMJ(KP,1)/AL(KP)
  111.             BMY=UMJ(KP,2)/AH(KP)
  112.             BM(KP)=BMX*BMX+BMY*BMY
  113.             BM(KP)=SQRT(BM(KP))+xpeti2
  114.             BPX=UPI(KP,1)/AL(KP)
  115.             BPY=UPI(KP,2)/AH(KP)
  116.             BP(KP)=BPX*BPX+BPY*BPY
  117.             BP(KP)=SQRT(BP(KP))+xpeti2
  118. 70074    CONTINUE
  119. C
  120. C     DECENTREMENT
  121. C
  122.          IF(IDCEN.EQ.1)THEN
  123. C CENTREE
  124.             DO 70081 K=KDEB,KFIN
  125.                KP=K-KDEB+1
  126.                TO1(KP)=0.D0
  127.                TO2(KP)=0.D0
  128.                CXT(KP)=0.D0
  129.                CYT(KP)=0.D0
  130.                CXY(KP)=0.D0
  131. 70081       CONTINUE
  132.          ELSEIF(IDCEN.EQ.2)THEN
  133. C SUPGDC
  134.             DO 7008 K=KDEB,KFIN
  135.                KP=K-KDEB+1
  136.                HMK=2.D0*UM(KP)/BM(KP)
  137.                XMB(KP)=HMK
  138.                ALFA=UM(KP)*HMK/COEFT(KP)/2.D0
  139.                AKSI=ALFA/3.D0
  140.                IF(ALFA.GT.3.D0)AKSI=1.D0
  141.                CCT=AKSI/BM(KP)
  142.                HPK=2.D0*UP(KP)/BP(KP)
  143.                ALFA=UP(KP)*HPK/COEFT(KP)/2.D0
  144.                AKSI=ALFA/3.D0
  145.                IF(ALFA.GT.3.D0)AKSI=1.D0
  146.                CCP=AKSI/BP(KP)
  147.                CPT=CCP-CCT
  148.                CC2=0.D0
  149.                IF(CPT.GE.0.D0)CC2=CPT
  150.                CC2=CC2*CC0
  151.                TO1(KP)=CCT
  152.                TO2(KP)=CC2
  153.                CXT(KP)=(UMJ(KP,1)*UMJ(KP,1)*CCT+UPI(KP,1)*UPI(KP,1)*CC2)
  154.                CYT(KP)=(UMJ(KP,2)*UMJ(KP,2)*CCT+UPI(KP,2)*UPI(KP,2)*CC2)
  155.                CXY(KP)=(UMJ(KP,1)*UMJ(KP,2)*CCT+UPI(KP,1)*UPI(KP,2)*CC2)
  156.  7008       CONTINUE
  157.          ELSEIF(IDCEN.EQ.3)THEN
  158. C SUPG
  159.             DO 7018 K=KDEB,KFIN
  160.                KP=K-KDEB+1
  161.                HMK=2.D0*UM(KP)/BM(KP)
  162.                XMB(KP)=HMK
  163.                ALFA=UM(KP)*HMK/COEFT(KP)/2.D0
  164.                AKSI=ALFA/3.D0
  165.                IF(ALFA.GT.3.D0)AKSI=1.D0
  166.                CCT=AKSI/BM(KP)
  167.                TO1(KP)=CCT
  168.                TO2(KP)=0.D0
  169.                CXT(KP)=UMJ(KP,1)*UMJ(KP,1)*CCT
  170.                CYT(KP)=UMJ(KP,2)*UMJ(KP,2)*CCT
  171.                CXY(KP)=UMJ(KP,1)*UMJ(KP,2)*CCT
  172.  7018       CONTINUE
  173.          ELSEIF(IDCEN.EQ.4)THEN
  174. C Tenseur Visqueux
  175.             DT19=DTM1*0.5D0
  176.             DO 7009 K=KDEB,KFIN
  177.                KP=K-KDEB+1
  178.                TO1(KP)=DT19
  179.                TO2(KP)=0.D0
  180.                CXT(KP)=UMJ(KP,1)*UMJ(KP,1)*DT19
  181.                CYT(KP)=UMJ(KP,2)*UMJ(KP,2)*DT19
  182.                CXY(KP)=UMJ(KP,1)*UMJ(KP,2)*DT19
  183.  7009       CONTINUE
  184.          ENDIF
  185. C
  186. C     CALCUL DU PAS DE TEMPS
  187. C
  188. C*IBMDIR* PREFER SCALAR
  189.          DO 7010 K=KDEB,KFIN
  190.             KP=K-KDEB+1
  191.             DT0=DT
  192.             DT1=2.D0/
  193.      &           (UMJ(KP,1)*UMJ(KP,1)/(CXT(KP)+COEFT(KP))
  194.      &           +UMJ(KP,2)*UMJ(KP,2)/(CYT(KP)+COEFT(KP)))
  195.             DT2=0.5D0/
  196.      &           ((CXT(KP)+COEFT(KP))*AL2(KP)
  197.      &           +(CYT(KP)+COEFT(KP))*AH2(KP))
  198.             IF(DT1.LT.DT)DT=DT1
  199.             IF(DT2.LT.DT)DT=DT2
  200.             IF(DT.NE.DT0) THEN
  201.                DTT1=DT1
  202.                DTT2=DT2
  203.                DIAEL=XMB(KP)
  204.                NUEL=K
  205.             END IF
  206.  7010    CONTINUE
  207.          DO 7011 I=1,NP
  208.             DO 70111 K=KDEB,KFIN
  209.                KP=K-KDEB+1
  210.                WTILDE(KP,I)=DRR(I,K)
  211.      &              +(TO1(KP)*(UMJ(KP,1)*HR(K,I,1)+UMJ(KP,2)*HR(K,I,2))
  212.      &              + TO2(KP)*(UPI(KP,1)*HR(K,I,1)+UPI(KP,2)*HR(K,I,2))
  213.      &              )*AIRE(KP)
  214. 70111       CONTINUE
  215.  7011    CONTINUE
  216.          DO 7012 K=KDEB,KFIN
  217.             KP=K-KDEB+1
  218.             CHGLD(KP)=COEFT(KP)/12.D0*QUA4*(AL2(KP)+AH2(KP))
  219.             CHGLP(KP)=(CXT(KP)*AL2(KP)+CYT(KP)*AH2(KP))/12.D0*QUA4
  220.  7012    CONTINUE
  221. C                 ***********
  222. C                 *   3 D   *
  223. C                 ***********
  224.       ELSEIF(IDIM.EQ.3)THEN
  225.          CUB8=CHGL
  226.          DO 8004 K=KDEB,KFIN
  227.             KP=K-KDEB+1
  228.             AL2(KP)=1.D0/AL(KP)/AL(KP)
  229.             AH2(KP)=1.D0/AH(KP)/AH(KP)
  230.             AP2(KP)=1.D0/AP(KP)/AP(KP)
  231.             XMB(KP)=(AL(KP)+AH(KP)+AP(KP))/3.D0
  232.             CFM(KP)=AL(KP)*AH(KP)/AP(KP)+AL(KP)*AP(KP)/AH(KP)+
  233.      &           AP(KP)*AH(KP)/AL(KP)
  234.  8004    CONTINUE
  235.          DO 8007 K=KDEB,KFIN
  236.             KP=K-KDEB+1
  237.             UMI(KP,1)=UMI(KP,1)/AIRE(KP)
  238.             UMI(KP,2)=UMI(KP,2)/AIRE(KP)
  239.             UMI(KP,3)=UMI(KP,3)/AIRE(KP)
  240.             UM(KP)=UMI(KP,1)*UMI(KP,1)+UMI(KP,2)*UMI(KP,2)
  241.      &           +UMI(KP,3)*UMI(KP,3)
  242.             UM(KP)=SQRT(UM(KP))
  243.             AX=GRAD(KP,1)*GRAD(KP,1)+GRAD(KP,2)*GRAD(KP,2)
  244.      &           +GRAD(KP,3)*GRAD(KP,3)
  245.             AX=SQRT(AX)+xpeti2
  246.             GRAD(KP,1)=GRAD(KP,1)/AX
  247.             GRAD(KP,2)=GRAD(KP,2)/AX
  248.             GRAD(KP,3)=GRAD(KP,3)/AX
  249.             UP(KP)=GRAD(KP,1)*UMI(KP,1)+GRAD(KP,2)*UMI(KP,2)
  250.      &           +GRAD(KP,3)*UMI(KP,3)
  251.             UPI(KP,1)=UP(KP)*GRAD(KP,1)
  252.             UPI(KP,2)=UP(KP)*GRAD(KP,2)
  253.             UPI(KP,3)=UP(KP)*GRAD(KP,3)
  254.  8007    CONTINUE
  255.          DO 80074 K=KDEB,KFIN
  256.             KP=K-KDEB+1
  257.             BMX=UMI(KP,1)/AL(KP)
  258.             BMY=UMI(KP,2)/AH(KP)
  259.             BMZ=UMI(KP,3)/AP(KP)
  260.             BM(KP)=BMX*BMX+BMY*BMY+BMZ*BMZ
  261.             BM(KP)=SQRT(BM(KP))+xpeti2
  262.             BPX=UPI(KP,1)/AL(KP)
  263.             BPY=UPI(KP,2)/AH(KP)
  264.             BPZ=UPI(KP,3)/AP(KP)
  265.             BP(KP)=BPX*BPX+BPY*BPY+BPZ*BPZ
  266.             BP(KP)=SQRT(BP(KP))+xpeti2
  267. 80074    CONTINUE
  268. C
  269. C     DECENTREMENT
  270. C
  271.          IF(IDCEN.EQ.1)THEN
  272. C CENTREE
  273.             DO 80081 K=KDEB,KFIN
  274.                KP=K-KDEB+1
  275.                TO1(KP)=0.D0
  276.                TO2(KP)=0.D0
  277.                CXT(KP)=0.D0
  278.                CYT(KP)=0.D0
  279.                CZT(KP)=0.D0
  280.                CXY(KP)=0.D0
  281.                CXZ(KP)=0.D0
  282.                CYZ(KP)=0.D0
  283. 80081       CONTINUE
  284.          ELSEIF(IDCEN.EQ.2)THEN
  285. C SUPGDC
  286.             DO 8008 K=KDEB,KFIN
  287.                KP=K-KDEB+1
  288.                HMK=2.D0*UM(KP)/BM(KP)
  289.                XMB(KP)=HMK
  290.                ALFA=UM(KP)*HMK/COEFT(KP)/2.D0
  291.                AKSI=ALFA/3.D0
  292.                IF(ALFA.GT.3.D0)AKSI=1.D0
  293.                CCT=AKSI/BM(KP)
  294.                HPK=2.D0*UP(KP)/BP(KP)
  295.                ALFA=UP(KP)*HPK/COEFT(KP)/2.D0
  296.                AKSI=ALFA/3.D0
  297.                IF(ALFA.GT.3.D0)AKSI=1.D0
  298.                CCP=AKSI/BP(KP)
  299.                CPT=CCP-CCT
  300.                CC2=0.D0
  301.                IF(CPT.GE.0.D0)CC2=CPT
  302.                CC2=CC2*CC0
  303.                TO1(KP)=CCT
  304.                TO2(KP)=CC2
  305.                CXT(KP)=(UMI(KP,1)*UMI(KP,1)*CCT+UPI(KP,1)*UPI(KP,1)*CC2)
  306.                CYT(KP)=(UMI(KP,2)*UMI(KP,2)*CCT+UPI(KP,2)*UPI(KP,2)*CC2)
  307.                CZT(KP)=(UMI(KP,3)*UMI(KP,3)*CCT+UPI(KP,3)*UPI(KP,3)*CC2)
  308.                CXY(KP)=(UMI(KP,1)*UMI(KP,2)*CCT+UPI(KP,1)*UPI(KP,2)*CC2)
  309.                CXZ(KP)=(UMI(KP,1)*UMI(KP,3)*CCT+UPI(KP,1)*UPI(KP,3)*CC2)
  310.                CYZ(KP)=(UMI(KP,2)*UMI(KP,3)*CCT+UPI(KP,2)*UPI(KP,3)*CC2)
  311.  8008       CONTINUE
  312.          ELSEIF(IDCEN.EQ.3)THEN
  313. C SUPG
  314.             DO 8018 K=KDEB,KFIN
  315.                KP=K-KDEB+1
  316.                HMK=2.D0*UM(KP)/BM(KP)
  317.                XMB(KP)=HMK
  318.                ALFA=UM(KP)*HMK/COEFT(KP)/2.D0
  319.                AKSI=ALFA/3.D0
  320.                IF(ALFA.GT.3.D0)AKSI=1.D0
  321.                CCT=AKSI/BM(KP)
  322.                TO1(KP)=CCT
  323.                TO2(KP)=0.D0
  324.                CXT(KP)=UMI(KP,1)*UMI(KP,1)*CCT
  325.                CYT(KP)=UMI(KP,2)*UMI(KP,2)*CCT
  326.                CZT(KP)=UMI(KP,3)*UMI(KP,3)*CCT
  327.                CXY(KP)=UMI(KP,1)*UMI(KP,2)*CCT
  328.                CXZ(KP)=UMI(KP,1)*UMI(KP,3)*CCT
  329.                CYZ(KP)=UMI(KP,2)*UMI(KP,3)*CCT
  330.  8018       CONTINUE
  331.          ELSEIF(IDCEN.EQ.4)THEN
  332. C Tenseur Visqueux
  333.             DT19=DTM1*0.5D0
  334.             DO 8009 K=KDEB,KFIN
  335.                KP=K-KDEB+1
  336.                TO1(KP)=DT19
  337.                TO2(KP)=0.D0
  338.                CXT(KP)=UMI(KP,1)*UMI(KP,1)*DT19
  339.                CYT(KP)=UMI(KP,2)*UMI(KP,2)*DT19
  340.                CZT(KP)=UMI(KP,3)*UMI(KP,3)*DT19
  341.                CXY(KP)=UMI(KP,1)*UMI(KP,2)*DT19
  342.                CXZ(KP)=UMI(KP,1)*UMI(KP,3)*DT19
  343.                CYZ(KP)=UMI(KP,2)*UMI(KP,3)*DT19
  344.  8009       CONTINUE
  345.          ENDIF
  346. C
  347. C     CALCUL DU PAS DE TEMPS
  348. C
  349. C*IBMDIR* PREFER SCALAR
  350.          DO 8010 K=KDEB,KFIN
  351.             KP=K-KDEB+1
  352.             DT0=DT
  353.             DT1=2.D0/
  354.      &           (UMI(KP,1)*UMI(KP,1)/(CXT(KP)+COEFT(KP))
  355.      &           +UMI(KP,2)*UMI(KP,2)/(CYT(KP)+COEFT(KP))
  356.      &           +UMI(KP,3)*UMI(KP,3)/(CZT(KP)+COEFT(KP)))
  357.             DT2=0.5D0/
  358.      &           ((CXT(KP)+COEFT(KP))*AL2(KP)
  359.      &           +(CYT(KP)+COEFT(KP))*AH2(KP)
  360.      &           +(CZT(KP)+COEFT(KP))*AP2(KP))
  361.             IF(DT1.LT.DT)DT=DT1
  362.             IF(DT2.LT.DT)DT=DT2
  363.             IF(DT.NE.DT0) THEN
  364.                DTT1=DT1
  365.                DTT2=DT2
  366.                DIAEL=XMB(KP)
  367.                NUEL=K
  368.             END IF
  369.  8010    CONTINUE
  370.          DO 8011 I=1,NP
  371.             DO 80111 K=KDEB,KFIN
  372.                KP=K-KDEB+1
  373.                WTILDE(KP,I)=DRR(I,K)
  374.      &              +(TO1(KP)*(UMI(KP,1)*HR(K,I,1)+UMI(KP,2)*HR(K,I,2)
  375.      &              +UMI(KP,3)*HR(K,I,3))
  376.      &              + TO2(KP)*(UPI(KP,1)*HR(K,I,1)+UPI(KP,2)*HR(K,I,2)
  377.      &              +UPI(KP,3)*HR(K,I,3))
  378.      &              )*AIRE(KP)
  379. 80111       CONTINUE
  380.  8011    CONTINUE
  381.          DO 8012 K=KDEB,KFIN
  382.             KP=K-KDEB+1
  383.             GEO1=CFM(KP)*CUB8
  384.             CHGLD(KP)=COEFT(KP)*GEO1
  385.             CHGLP(KP)=(CXT(KP)+CYT(KP)+CZT(KP))/3.D0*GEO1
  386.  8012    CONTINUE
  387.       ENDIF
  388.       RETURN
  389.  1002 FORMAT(10(1X,1PE11.4))
  390.       END
  391.  
  392.  
  393.  
  394.  
  395.  
  396.  
  397.  
  398.  
  399.  
  400.  
  401.  

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