Télécharger la plaquette Cast3M


fr en

Télécharger zpsi.eso

Retour à la liste

Numérotation des lignes : 

zpsi

  1. C ZPSI      SOURCE    CHAT      05/01/13    04:24:52     5004
  2.       SUBROUTINE ZPSI
  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.      & (HR,RPG,DRR,LE,NEL,K0,NPTD,IES,NP,IAXI,
  12.      & IPADL,IKOMP,IKAS,
  13.      & ALFE,IND1,UN,INDU,NPTU,IPADU,
  14.      & TN,QE,IKS,
  15.      & HRN,G,ALT,SGT,
  16.      & VOLU,COTE,NELZ,ZTE,
  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. -INC CCVQUA4
  37. -INC CCREEL
  38. -INC SMCOORD
  39.  
  40. -INC PPARAM
  41. -INC CCOPTIO
  42.  
  43. C
  44. C Longueur des registres vectoriels de la machine cible
  45. C On prend 64 pour ne pas augmenter la taille des tableaux
  46. C nécessaires à la vectorisation.
  47. C
  48.       PARAMETER(LRV=64)
  49.  
  50.       DIMENSION UN(NPTU,IES),HRN(NPTD),TN(NPTD)
  51.       DIMENSION COTE(NELZ,IES),VOLU(NELZ),QE(*)
  52.       DIMENSION ALFE(*),ALT(*)
  53.  
  54.       DIMENSION IPADL(1),LE(NP,1),IPADU(*)
  55.       DIMENSION HR(NEL,NP,IES),RPG(1),DRR(NP,NEL)
  56.  
  57.       DIMENSION BF(9,9)
  58.       DIMENSION QGGT(8,8),Q1(8,8),Q2(8,8),Q3(8,8)
  59.  
  60.       DIMENSION COEF(LRV),AIRE(LRV)
  61.       DIMENSION AL2(LRV),AH2(LRV),AP2(LRV)
  62.       DIMENSION AL(LRV),AH(LRV),AP(LRV)
  63.       DIMENSION ALFT(LRV),QT(LRV)
  64.       DIMENSION XMB(LRV),XMH(LRV)
  65.       DIMENSION CFM(LRV)
  66.       DIMENSION CF1(LRV),CF2(LRV),CF3(LRV)
  67.       DIMENSION DR(LRV,9)
  68.       DIMENSION UM(LRV),UP(LRV)
  69.       DIMENSION UIX(LRV,9),UIY(LRV,9),UIZ(LRV,9)
  70.       DIMENSION TETAC(LRV,9),TETAD(LRV,9)
  71.       DIMENSION UMI(LRV,3)
  72.       DIMENSION SBF(LRV,9)
  73.       DIMENSION GRADT(LRV,3)
  74.  
  75. C?      DIMENSION ZTE(LRVH)
  76.       DIMENSION ZTE(NPTU)
  77.  
  78.       INTEGER   U(3),D(3),SGN
  79.       REAL*8    G(1),ZVGG(4),x(4),y(4),KMAX
  80.       REAL*8    KKK(4),PhiSource(LRV),n(3,4)
  81.       REAL*8    PhiN3,PhiN4,Phi3NNQ,Phi4NNQ
  82.       REAL*8    MINI2,MINI4,MAXI2,MAXI4
  83.  
  84.       SAVE      IPAS,QGGT,Q1,Q2,Q3
  85.  
  86.       DATA      CD/1.D0/
  87.       DATA      IPAS/0/
  88. C*********************************************************************
  89. C
  90. C           INITIALISATIONS DIVERSES
  91. C
  92.  
  93.       NK=K0
  94. C                           ********
  95. C                           *  2D  *
  96. C                           ********
  97.  
  98.       IF (IES.EQ.3) GOTO 10
  99.  
  100.       IAX1=0
  101.       IAX2=0
  102.       IF (IAXI.EQ.1) IAX2=1
  103.       IF (IAXI.EQ.2) IAX1=1
  104.  
  105. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  106.  
  107. C                           DIFFERENCES TRIANGLE / QUADRANGLE
  108.       QUA4=0.D0
  109.       IF (NP.EQ.4) QUA4=1.D0
  110. C
  111. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  112.  
  113. C
  114. C Calcul du nombre de paquets de LRV éléments
  115. C
  116.       NNN=MOD(NEL,LRV)
  117.       IF(NNN.EQ.0) NPACK=NEL/LRV
  118.       IF(NNN.NE.0) NPACK=1+(NEL-NNN)/LRV
  119.       KPACKD=1
  120.       KPACKF=NPACK
  121.  
  122. ******* BOUCLE SUR LES PAQUETS DE LRV ELEMENTS **********
  123.  
  124.       DO 7001 KPACK=KPACKD,KPACKF
  125.  
  126. C ======= A L'INTERIEUR DE CHAQUE PAQUET DE LRV ELEMENTS =======
  127.  
  128. C 1. Calcul des limites du paquet courant.
  129.  
  130.       KDEB=1+(KPACK-1)*LRV
  131.       KFIN=MIN(NEL,KDEB+LRV-1)
  132.       DO 7002 K=KDEB,KFIN
  133.       KP=K-KDEB+1
  134.       NK=K+K0
  135.       NK1=(1-IND1)*(NK-1)+1
  136.       ALFT(KP)=ALFE(NK1)+XPETIT
  137.       AIRE(KP)=VOLU(NK)
  138.       AL(KP)=COTE(NK,1)
  139.       AH(KP)=COTE(NK,2)
  140.       AL2(KP)=1.D0/AL(KP)/AL(KP)
  141.       AH2(KP)=1.D0/AH(KP)/AH(KP)
  142. 7002  CONTINUE
  143.  
  144.       IF((IKOMP.EQ.0.AND.IKAS.EQ.5).OR.
  145.      &(IKOMP.EQ.1.AND.IKAS.EQ.6))THEN
  146.       DO 7003 K=KDEB,KFIN
  147.       KP=K-KDEB+1
  148.       NK=K+K0
  149.       ALFT(KP)=ALFT(KP)+ALT(NK)/SGT
  150. 7003  CONTINUE
  151.       ENDIF
  152.  
  153. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  154.  
  155. C Initialisation des UMI avant accumulation
  156.       DO 7005 K=KDEB,KFIN
  157.       KP=K-KDEB+1
  158.       UMI(KP,1)=XPETIT
  159.       UMI(KP,2)=XPETIT
  160. 7005  CONTINUE
  161.  
  162.       DO 7006 I=1,NP
  163.       DO 7006 K=KDEB,KFIN
  164.       KP=K-KDEB+1
  165.       NF=IPADL(LE(I,K))
  166.       NFU=IPADU(LE(I,K))
  167.       NFU=(1-INDU)*(NFU-1)+1
  168.       DR(KP,I)=DRR(I,K)
  169.       UMI(KP,1)=UMI(KP,1)+UN(NFU,1)*DR(KP,I)
  170.       UMI(KP,2)=UMI(KP,2)+UN(NFU,2)*DR(KP,I)
  171.  
  172.       TETAC(KP,I)=HRN(NF)
  173.       TETAD(KP,I)=TN(NF)
  174. 7006  CONTINUE
  175.  
  176. *     Calcul du FLUX PhiSource du terme source.
  177.  
  178.       IF (IKOMP.EQ.0) THEN
  179.  
  180.          DO 70021 K=KDEB,KFIN
  181.          KP=K-KDEB+1
  182.          NK=K+K0
  183.          NKS=(1-IKS)*(NK-1)+1
  184.          QT(KP)=QE(NKS)
  185. 70021    CONTINUE
  186.  
  187.          DO 70062 K=KDEB,KFIN
  188.          KP=K-KDEB+1
  189.          PhiSource(KP)=-QT(KP)*AIRE(KP)
  190. 70062    CONTINUE
  191.  
  192.       ELSEIF (IKOMP.EQ.1) THEN
  193.  
  194.          DO 70023 K=KDEB,KFIN
  195.          KP=K-KDEB+1
  196.          NK=K+K0
  197.          NKS=(1-IKS)*(NK-1)+1
  198.          QT(KP)=QE(NKS)
  199. 70023    CONTINUE
  200.  
  201.          DO 70064 K=KDEB,KFIN
  202.          KP=K-KDEB+1
  203.          PhiSource(KP)=-QT(KP)*AIRE(KP)
  204. 70064    CONTINUE
  205.  
  206.       ENDIF
  207.  
  208.  
  209.       DO 7007 K=KDEB,KFIN
  210.       KP=K-KDEB+1
  211.       UMI(KP,1)=UMI(KP,1)/AIRE(KP)
  212.       UMI(KP,2)=UMI(KP,2)/AIRE(KP)
  213.       UM(KP)=UMI(KP,1)*UMI(KP,1)+UMI(KP,2)*UMI(KP,2)
  214.       UM(KP)=SQRT(UM(KP)) + XPETIT
  215.  
  216. 7007  CONTINUE
  217.  
  218. **************************************************
  219.  
  220.       CALL INITD(ZTE,NPTU,XPETIT)
  221.  
  222.       IF (IKOMP.EQ.0) THEN
  223.  
  224.       DO 7008 K=KDEB,KFIN
  225.       KP=K-KDEB+1
  226.  
  227.       DO 23 I=1,NP
  228. 23    ZVGG(I) =0.D0
  229.  
  230. **************************************************
  231.  
  232. * NP=3 correspond à PSI sur triangles
  233. * NP=4 correspond à PSI sur quadrangles ,cad NNQ
  234.  
  235. * Calcul des Ki:**********************************
  236.  
  237.       IF (NP.EQ.3) THEN
  238.  
  239.       DO 1 I=1,NP
  240. 1     KKK(I)=AIRE(KP)*(UMI(KP,1)*HR(K,I,1)+UMI(KP,2)*HR(K,I,2))
  241.  
  242.       ELSEIF (NP.EQ.4) THEN
  243.  
  244. * Calcul des normales:
  245.  
  246.       DO 36 I=1,NP
  247.  
  248.       x(I)=XCOOR((LE(I,K)-1)*(IES+1)+1)
  249. 36    y(I)=XCOOR((LE(I,K)-1)*(IES+1)+2)
  250.  
  251.       DO 39 I=1,NP
  252.  
  253.       IF (I.EQ.4) THEN
  254.       J1=1
  255.       ELSE
  256.       J1=I+1
  257.       ENDIF
  258.  
  259.       n(1,I)= (y(I)-y(J1))
  260.       n(2,I)= (x(J1)-x(I))
  261.  
  262. 39    CONTINUE
  263.  
  264. * Calcul des Ki :
  265.  
  266.       DO 17 I=1,NP
  267. 17    KKK(I)=0.5D0*(UMI(KP,1)*n(1,I)+UMI(KP,2)*n(2,I))
  268.  
  269.       ENDIF
  270.  
  271.  
  272. ****************************************************
  273. * Calcul du DT:
  274. ****************************************************
  275.  
  276.       KMAX=KKK(1)
  277.       DO 29 I=2,NP
  278.       IF (KKK(I).GT.KMAX) KMAX=KKK(I)
  279. 29    CONTINUE
  280.  
  281.       DO 27 J=1,NP
  282.       IP=IPADU(LE(J,K))
  283.  
  284.       ZTE(IP) = ZTE(IP) + KMAX/(2.D0*AIRE(KP)+XPETIT)
  285. 27    CONTINUE
  286.  
  287. ****************************************************
  288.  
  289.  
  290.       IF (NP.EQ.4) THEN
  291.  
  292. **************************************************
  293. *     Schéma NNQ (Pour des Quadrangles) .
  294. **************************************************
  295.  
  296. * Tests:
  297. ********
  298.       Nd=0
  299.  
  300.       rrt =KKK(2)+KKK(3)
  301.  
  302.       IF (rrt.GT.0.D0) THEN
  303.       Nd=Nd+1
  304.       D(Nd)=1
  305.       ENDIF
  306.  
  307.       rrt =KKK(3)+KKK(4)
  308.  
  309.       IF (rrt.GT.0.D0) THEN
  310.       Nd=Nd+1
  311.       D(Nd)=2
  312.       ENDIF
  313.  
  314.       rrt =KKK(1)+KKK(4)
  315.  
  316.       IF (rrt.GT.0.D0) THEN
  317.       Nd=Nd+1
  318.       D(Nd)=3
  319.       ENDIF
  320.  
  321.       rrt =KKK(1)+KKK(2)
  322.  
  323.       IF (rrt.GT.0.D0) THEN
  324.       Nd=Nd+1
  325.       D(Nd)=4
  326.       ENDIF
  327.  
  328.       IF (Nd.EQ.2) THEN
  329.       IF ((D(1)+1).NE.D(2)) THEN
  330.       jj = D(1)
  331.       D(1) = D(2)
  332.       D(2) = jj
  333.       ENDIF
  334.       ENDIF
  335.  
  336.       Nu=0
  337.       DO 66 I=1,NP
  338.       IF (Nd.EQ.2) THEN
  339.       IF ((I.NE.D(1)).AND.(I.NE.D(2))) THEN
  340.       Nu=Nu+1
  341.       U(Nu)=I
  342.       ENDIF
  343.       ELSEIF (Nd.EQ.1) THEN
  344.       IF (I.NE.D(1)) THEN
  345.       Nu=Nu+1
  346.       U(Nu)=I
  347.       ENDIF
  348.       ENDIF
  349. 66    CONTINUE
  350.  
  351.       IF (Nu+Nd.NE.4) THEN
  352.       Print *,'Nd=,Nu=',Nd,Nu
  353.       ENDIF
  354.  
  355.       IF (Nu.EQ.2) THEN
  356.       IF ((U(1)+1).NE.U(2)) THEN
  357.       jj=U(1)
  358.       U(1)=U(2)
  359.       U(2)=jj
  360.       ENDIF
  361.       ENDIF
  362.  
  363.       N1 =U(1)
  364.       N2 =U(2)
  365.       N3 =D(1)
  366.       N4 =D(2)
  367.  
  368.       IF (Nd.EQ.2) THEN
  369.  
  370.       MINI2=0.D0
  371.       IF (0.D0.GE.KKK(N2))      MINI2=KKK(N2)
  372.  
  373.       MINI4=0.D0
  374.       IF (0.D0.GE.KKK(N4))      MINI4=KKK(N4)
  375.  
  376.       MAXI2=0.D0
  377.       IF (0.D0.LE.KKK(N2))      MAXI2=KKK(N2)
  378.  
  379.       MAXI4=0.D0
  380.       IF (0.D0.LE.KKK(N4))      MAXI4=KKK(N4)
  381.  
  382.       PhiN3 =( KKK(N1) + MINI2 + MINI4 )
  383.      *  *(TETAC(KP,N3)-TETAC(KP,N2))
  384.      *  +( MAXI4-MINI2 )*(TETAC(KP,N3)-TETAC(KP,N1))
  385.  
  386.       PhiN4 =( KKK(N1) + MINI2 + MINI4 )
  387.      * *(TETAC(KP,N4)-TETAC(KP,N1))
  388.      * +( MAXI2-MINI4 )*(TETAC(KP,N4)-TETAC(KP,N2))
  389.  
  390.       ZVGG(N3)=PhiN3
  391.       ZVGG(N4)=PhiN4
  392.  
  393.       r= -ZVGG(N3) / (ZVGG(N4)+XPETIT)
  394.       IF (r.GT.0.D0) THEN
  395.       IF (r.GT.1.D0) THEN
  396.       ZVGG(N3)=ZVGG(N3)+ZVGG(N4)
  397.       ZVGG(N4)=0.D0
  398.       ELSE
  399.       ZVGG(N4)=ZVGG(N3)+ZVGG(N4)
  400.       ZVGG(N3)=0.D0
  401.       ENDIF
  402.       ENDIF
  403.  
  404.       ELSEIF (ND.EQ.1) THEN
  405.  
  406.       PhiQ =0.5D0*(UMI(KP,1)*(n(1,1)+n(1,4)) +
  407.      * UMI(KP,2)*(n(2,1)+n(2,4)))*(TETAC(KP,3)-TETAC(KP,1))+
  408.      * 0.5D0*(UMI(KP,1)*(n(1,1)+n(1,2)) +
  409.      * UMI(KP,2)*(n(2,1)+n(2,2)))*(TETAC(KP,4)-TETAC(KP,2))
  410.  
  411.       ZVGG(N3) = PhiQ
  412.       ENDIF
  413.  
  414.       ELSEIF (NP.EQ.3) THEN
  415.  
  416. **************************************************
  417. *     Schéma PSI (Min-Mod) Pour des Triangles :
  418. **************************************************
  419.  
  420. *  U Noeuds amonts
  421. *  D Noeuds avals .
  422. **********************
  423.  
  424.       Nd=0
  425.       Nu=0
  426.  
  427.       DO 13 I=1,NP
  428.  
  429.       IF (KKK(I).GT.0.D0) THEN
  430.       Nd=Nd+1
  431.       D(Nd)=I
  432.       ELSE
  433.       Nu=Nu+1
  434.       U(Nu)=I
  435.       ENDIF
  436.  
  437. 13    CONTINUE
  438.  
  439.       IF (Nd.EQ.1) THEN
  440.  
  441. ***************************
  442. * 1 Target :
  443. ***************************
  444.  
  445.       N1 = D(1)
  446.       N2 = U(1)
  447.       N3 = U(2)
  448.  
  449.       ZVGG(N1) = PhiSource(KP) + KKK(1)*TETAC(KP,1)+
  450.      * KKK(2)*TETAC(KP,2) + KKK(3)*TETAC(KP,3)
  451.       ZVGG(N2) = 0.D0
  452.       ZVGG(N3) = 0.D0
  453.  
  454.       ELSEIF (Nd.EQ.2) THEN
  455.  
  456. ***************************
  457. * 2 Targets :
  458. ***************************
  459.  
  460.       N1 = D(1)
  461.       N2 = D(2)
  462.       N3 = U(1)
  463.  
  464.       ZVGG(N1) = KKK(N1)*(TETAC(KP,N1)-TETAC(KP,N3))-PhiSource(KP)*
  465.      * KKK(N1)/(KKK(N3)+XPETIT)
  466.  
  467.       ZVGG(N2) = KKK(N2)*(TETAC(KP,N2)-TETAC(KP,N3))-PhiSource(KP)*
  468.      * KKK(N2)/(KKK(N3)+XPETIT)
  469.  
  470.       ZVGG(N3) = 0.D0
  471.  
  472.       r= -ZVGG(N1) / (ZVGG(N2)+XPETIT)
  473.       IF (r.GT.0.D0) THEN
  474.       IF (r.GT.1.D0) THEN
  475.       ZVGG(N1)=ZVGG(N1)+ZVGG(N2)
  476.       ZVGG(N2)=0.D0
  477.       ELSE
  478.       ZVGG(N2)=ZVGG(N1)+ZVGG(N2)
  479.       ZVGG(N1)=0.D0
  480.       ENDIF
  481.       ENDIF
  482.  
  483.       ENDIF
  484.  
  485.       ENDIF
  486.  
  487.       DO 7 I=1,NP
  488.       SBF(KP,I)=ZVGG(I)
  489. 7     CONTINUE
  490.  
  491. 7008  CONTINUE
  492. **************************************************
  493.  
  494.       ELSEIF (IKOMP.EQ.1) THEN
  495.  
  496. * Prog à compléter
  497.       RETURN
  498.  
  499.       ENDIF
  500.  
  501. * Pas de temps DT :
  502.  
  503.       DT=100.
  504.  
  505.       DO 32 K1=KDEB,KFIN
  506.       DO 32 J=1,NP
  507.       IP=IPADU(LE(J,K1))
  508.       DT1=1.D0 /( ZTE(IP)+XPETIT)
  509.  
  510. c     IF (Abs(DT1).GE.1000) THEN
  511. c     Print *,'DT1 trop grand =',DT1
  512. c     ENDIF
  513.  
  514.       IF (DT1.LE.DT) DT=DT1
  515. 32    CONTINUE
  516.  
  517.  
  518.       DTT1=DT
  519.       DTT2=DT
  520.       DTT3=0.D0
  521.       DIAEL=0.D0
  522.       NUEL=0
  523.  
  524. *      PRINT *, 'PAS DE TEMPS = ', DT
  525.  
  526.  
  527. ****************************************************
  528. * Contribution du terme diffusif:
  529.  
  530.       DO 8 I=1,NP
  531.       DO 8 J=1,NP
  532.       DO 8 K=KDEB,KFIN
  533.       KP=K-KDEB+1
  534.  
  535.       ZVGT=AIRE(KP)*(
  536.      &     HR(K,I,1)*HR(K,J,1)*ALFT(KP)
  537.      &+    HR(K,I,2)*HR(K,J,2)*ALFT(KP)
  538.      &+    (ALFT(KP)*AL2(KP)+ALFT(KP)*AH2(KP))/12.D0
  539.      &     *VGGT(J,I)*QUA4  )
  540.  
  541.  
  542.       SBF(KP,I)=SBF(KP,I)+ TETAD(KP,J)*ZVGT
  543.  
  544. 8     CONTINUE
  545. ****************************************************
  546.  
  547. C     Fin de l'accumulation dans SBF.
  548. C     On ajoute ces incréments G.
  549.  
  550.  
  551.       DO 7017 I=1,NP
  552.       DO 7017 K=KDEB,KFIN
  553.       KP=K-KDEB+1
  554.       NF=IPADL(LE(I,K))
  555.  
  556.       G(NF) = G(NF)+SBF(KP,I)
  557.  
  558. 7017  CONTINUE
  559.  
  560.  
  561. 7001  CONTINUE
  562.       IPAS=1
  563.  
  564.       RETURN
  565.  
  566.  
  567. C                           ********
  568. C                           *  3D  *
  569. C                           ********
  570.  
  571.  10   CONTINUE
  572.  
  573.       IF (IPAS.EQ.0) CALL CALHRH(QGGT,Q1,Q2,Q3,IES)
  574.       CUB8 = 0.D0
  575.       IF (NP.EQ.8) CUB8=1.D0
  576.  
  577. C
  578. C Calcul du nombre de paquets de LRV éléments
  579. C
  580.       NNN=MOD(NEL,LRV)
  581.       IF(NNN.EQ.0) NPACK=NEL/LRV
  582.       IF(NNN.NE.0) NPACK=1+(NEL-NNN)/LRV
  583.       KPACKD=1
  584.       KPACKF=NPACK
  585. C
  586. C ******* BOUCLE SUR LES PAQUETS DE LRV ELEMENTS **********
  587. C
  588.       DO 8001 KPACK=KPACKD,KPACKF
  589. C
  590. C ======= A L'INTERIEUR DE CHAQUE PAQUET DE LRV ELEMENTS =======
  591. C
  592. C 1. Calcul des limites du paquet courant.
  593.       KDEB=1+(KPACK-1)*LRV
  594.       KFIN=MIN(NEL,KDEB+LRV-1)
  595. C
  596.       DO 8002 K=KDEB,KFIN
  597.       KP=K-KDEB+1
  598.       NK=K+K0
  599.       NK1=(1-IND1)*(NK-1)+1
  600.       ALFT(KP)=ALFE(NK1)+XPETIT
  601.       AIRE(KP)=VOLU(NK)
  602.       AL(KP)=COTE(NK,1)
  603.       AH(KP)=COTE(NK,2)
  604.       AP(KP)=COTE(NK,3)
  605.  
  606.       CFM(KP)=AL(KP)*AH(KP)/AP(KP)+AL(KP)*AP(KP)/AH(KP)+
  607.      &        AP(KP)*AH(KP)/AL(KP)
  608. C     CF1(KP)=AL(KP)*AH(KP)/AP(KP)
  609. C     CF2(KP)=AL(KP)*AP(KP)/AH(KP)
  610. C     CF3(KP)=AP(KP)*AH(KP)/AL(KP)
  611.       XMH(KP)=(AL(KP)+AH(KP)+AP(KP))/3.D0
  612.       AL2(KP)=1.D0/AL(KP)/AL(KP)
  613.       AH2(KP)=1.D0/AH(KP)/AH(KP)
  614.       AP2(KP)=1.D0/AP(KP)/AP(KP)
  615.  
  616.  8002 CONTINUE
  617.       IF((IKOMP.EQ.0.AND.IKAS.EQ.5).OR.
  618.      &(IKOMP.EQ.1.AND.IKAS.EQ.6))THEN
  619.       DO 8003 K=KDEB,KFIN
  620.       KP=K-KDEB+1
  621.       NK=K+K0
  622.       ALFT(KP)=ALFT(KP)+ALT(NK)/SGT
  623.  8003 CONTINUE
  624.       ENDIF
  625.  
  626. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  627.  
  628. C Initialisation des UMI avant accumulation
  629.       DO 8005 K=KDEB,KFIN
  630.       KP=K-KDEB+1
  631.       UMI(KP,1)=XPETIT
  632.       UMI(KP,2)=XPETIT
  633.       UMI(KP,3)=XPETIT
  634.  8005 CONTINUE
  635.  
  636.       DO 8006 I=1,NP
  637.       DO 8006 K=KDEB,KFIN
  638.       KP=K-KDEB+1
  639.       NF=IPADL(LE(I,K))
  640.       NFU=IPADU(LE(I,K))
  641.       NFU=(1-INDU)*(NFU-1)+1
  642.       DR(KP,I)=DRR(I,K)
  643.       UMI(KP,1)=UMI(KP,1)+UN(NFU,1)*DR(KP,I)
  644.       UMI(KP,2)=UMI(KP,2)+UN(NFU,2)*DR(KP,I)
  645.       UMI(KP,3)=UMI(KP,3)+UN(NFU,3)*DR(KP,I)
  646.  
  647.       TETAC(KP,I)=HRN(NF)
  648.       TETAD(KP,I)=TN(NF)
  649. 8006  CONTINUE
  650.  
  651. C
  652. C Initialisation de la variable d'accumulation SBF au terme source
  653. C
  654. C       write(6,*)' IKomp,ikas=',IKomp,ikas
  655. C       write(6,*)' IKS,IND1,INDU=',IKS,IND1,INDU
  656.  
  657.       IF(IKOMP.EQ.0)THEN
  658.  
  659.          DO 80021 K=KDEB,KFIN
  660.          KP=K-KDEB+1
  661.          NK=K+K0
  662.          NKS=(1-IKS)*(NK-1)+1
  663.          QT(KP)=QE(NKS)
  664. 80021    CONTINUE
  665.  
  666.          DO 80062 I=1,NP
  667.          DO 80062 K=KDEB,KFIN
  668.          KP=K-KDEB+1
  669.          SBF(KP,I)=-QT(KP)*DR(KP,I)
  670. 80062    CONTINUE
  671.  
  672.       ELSEIF(IKOMP.EQ.1)THEN
  673.  
  674.          DO 80023 K=KDEB,KFIN
  675.          KP=K-KDEB+1
  676.          NK=K+K0
  677.          NKS=(1-IKS)*(NK-1)+1
  678.          QT(KP)=QE(NKS)
  679. 80023    CONTINUE
  680.  
  681.          DO 80064 I=1,NP
  682.          DO 80064 K=KDEB,KFIN
  683.          KP=K-KDEB+1
  684.          SBF(KP,I)=-QT(KP)*DR(KP,I)
  685. 80064    CONTINUE
  686.  
  687.       ENDIF
  688.  
  689.       DO 8007 K=KDEB,KFIN
  690.       KP=K-KDEB+1
  691.       UMI(KP,1)=UMI(KP,1)/AIRE(KP)
  692.       UMI(KP,2)=UMI(KP,2)/AIRE(KP)
  693.       UMI(KP,3)=UMI(KP,3)/AIRE(KP)
  694.       UM(KP)=UMI(KP,1)*UMI(KP,1)+UMI(KP,2)*UMI(KP,2)
  695.      &      +UMI(KP,3)*UMI(KP,3)
  696.       UM(KP)=SQRT(UM(KP))
  697. 8007  CONTINUE
  698.  
  699.  
  700. ********************************************************
  701. *          Debut de PSI:
  702.  
  703.  
  704.  
  705.  
  706.  
  707.  
  708.  
  709.  
  710. ********************************************************
  711.  
  712.  
  713. C Le coeur du calcul ...
  714.  
  715. *      IF(IKOMP.EQ.0)THEN
  716.  
  717. C      DO 8014 I=1,NP
  718. C      DO 8014 K=KDEB,KFIN
  719. C      KP=K-KDEB+1
  720.  
  721. *      ELSEIF(IKOMP.EQ.1)THEN
  722. *      RETURN
  723.  
  724. *      ENDIF
  725.  
  726. C
  727. C     Fin de l'accumulation dans SBF.
  728. C     On ajoute ces incréments G.
  729. C
  730.       DO 8017 I=1,NP
  731.       DO 8017 K=KDEB,KFIN
  732.  
  733.       KP=K-KDEB+1
  734.       NF=IPADL(LE(I,K))
  735.       G(NF) = G(NF)+SBF(KP,I)
  736. 8017  CONTINUE
  737.  
  738. 8001  CONTINUE
  739.  
  740.  
  741.       IPAS=1
  742.       RETURN
  743. 1002  FORMAT(10(1X,1PE11.4))
  744.  
  745.       END
  746.  
  747.  
  748.  
  749.  
  750.  
  751.  

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