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    CB215821  26/08/24    21:19:05     12622          
  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 80065 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. 80065 CONTINUE
  176.  
  177. *     Calcul du FLUX PhiSource du terme source.
  178.  
  179.       IF (IKOMP.EQ.0) THEN
  180.  
  181.          DO 70021 K=KDEB,KFIN
  182.          KP=K-KDEB+1
  183.          NK=K+K0
  184.          NKS=(1-IKS)*(NK-1)+1
  185.          QT(KP)=QE(NKS)
  186. 70021    CONTINUE
  187.  
  188.          DO 70062 K=KDEB,KFIN
  189.          KP=K-KDEB+1
  190.          PhiSource(KP)=-QT(KP)*AIRE(KP)
  191. 70062    CONTINUE
  192.  
  193.       ELSEIF (IKOMP.EQ.1) THEN
  194.  
  195.          DO 70023 K=KDEB,KFIN
  196.          KP=K-KDEB+1
  197.          NK=K+K0
  198.          NKS=(1-IKS)*(NK-1)+1
  199.          QT(KP)=QE(NKS)
  200. 70023    CONTINUE
  201.  
  202.          DO 70064 K=KDEB,KFIN
  203.          KP=K-KDEB+1
  204.          PhiSource(KP)=-QT(KP)*AIRE(KP)
  205. 70064    CONTINUE
  206.  
  207.       ENDIF
  208.  
  209.  
  210.       DO 7007 K=KDEB,KFIN
  211.       KP=K-KDEB+1
  212.       UMI(KP,1)=UMI(KP,1)/AIRE(KP)
  213.       UMI(KP,2)=UMI(KP,2)/AIRE(KP)
  214.       UM(KP)=UMI(KP,1)*UMI(KP,1)+UMI(KP,2)*UMI(KP,2)
  215.       UM(KP)=SQRT(UM(KP)) + XPETIT
  216.  
  217. 7007  CONTINUE
  218.  
  219. **************************************************
  220.  
  221.       CALL INITD(ZTE,NPTU,XPETIT)
  222.  
  223.       IF (IKOMP.EQ.0) THEN
  224.  
  225.       DO 7008 K=KDEB,KFIN
  226.       KP=K-KDEB+1
  227.  
  228.       DO 23 I=1,NP
  229.       ZVGG(I) =0.D0
  230.    23 CONTINUE
  231.  
  232. **************************************************
  233.  
  234. * NP=3 correspond à PSI sur triangles
  235. * NP=4 correspond à PSI sur quadrangles ,cad NNQ
  236.  
  237. * Calcul des Ki:**********************************
  238.  
  239.       IF (NP.EQ.3) THEN
  240.  
  241.       DO 1 I=1,NP
  242.       KKK(I)=AIRE(KP)*(UMI(KP,1)*HR(K,I,1)+UMI(KP,2)*HR(K,I,2))
  243.     1 CONTINUE
  244.  
  245.       ELSEIF (NP.EQ.4) THEN
  246.  
  247. * Calcul des normales:
  248.  
  249.       DO 36 I=1,NP
  250.  
  251.       x(I)=XCOOR((LE(I,K)-1)*(IES+1)+1)
  252.       y(I)=XCOOR((LE(I,K)-1)*(IES+1)+2)
  253.    36 CONTINUE
  254.  
  255.       DO 39 I=1,NP
  256.  
  257.       IF (I.EQ.4) THEN
  258.       J1=1
  259.       ELSE
  260.       J1=I+1
  261.       ENDIF
  262.  
  263.       n(1,I)= (y(I)-y(J1))
  264.       n(2,I)= (x(J1)-x(I))
  265.  
  266. 39    CONTINUE
  267.  
  268. * Calcul des Ki :
  269.  
  270.       DO 17 I=1,NP
  271.       KKK(I)=0.5D0*(UMI(KP,1)*n(1,I)+UMI(KP,2)*n(2,I))
  272.    17 CONTINUE
  273.  
  274.       ENDIF
  275.  
  276.  
  277. ****************************************************
  278. * Calcul du DT:
  279. ****************************************************
  280.  
  281.       KMAX=KKK(1)
  282.       DO 29 I=2,NP
  283.       IF (KKK(I).GT.KMAX) KMAX=KKK(I)
  284. 29    CONTINUE
  285.  
  286.       DO 27 J=1,NP
  287.       IP=IPADU(LE(J,K))
  288.  
  289.       ZTE(IP) = ZTE(IP) + KMAX/(2.D0*AIRE(KP)+XPETIT)
  290. 27    CONTINUE
  291.  
  292. ****************************************************
  293.  
  294.  
  295.       IF (NP.EQ.4) THEN
  296.  
  297. **************************************************
  298. *     Schéma NNQ (Pour des Quadrangles) .
  299. **************************************************
  300.  
  301. * Tests:
  302. ********
  303.       Nd=0
  304.  
  305.       rrt =KKK(2)+KKK(3)
  306.  
  307.       IF (rrt.GT.0.D0) THEN
  308.       Nd=Nd+1
  309.       D(Nd)=1
  310.       ENDIF
  311.  
  312.       rrt =KKK(3)+KKK(4)
  313.  
  314.       IF (rrt.GT.0.D0) THEN
  315.       Nd=Nd+1
  316.       D(Nd)=2
  317.       ENDIF
  318.  
  319.       rrt =KKK(1)+KKK(4)
  320.  
  321.       IF (rrt.GT.0.D0) THEN
  322.       Nd=Nd+1
  323.       D(Nd)=3
  324.       ENDIF
  325.  
  326.       rrt =KKK(1)+KKK(2)
  327.  
  328.       IF (rrt.GT.0.D0) THEN
  329.       Nd=Nd+1
  330.       D(Nd)=4
  331.       ENDIF
  332.  
  333.       IF (Nd.EQ.2) THEN
  334.       IF ((D(1)+1).NE.D(2)) THEN
  335.       jj = D(1)
  336.       D(1) = D(2)
  337.       D(2) = jj
  338.       ENDIF
  339.       ENDIF
  340.  
  341.       Nu=0
  342.       DO 66 I=1,NP
  343.       IF (Nd.EQ.2) THEN
  344.       IF ((I.NE.D(1)).AND.(I.NE.D(2))) THEN
  345.       Nu=Nu+1
  346.       U(Nu)=I
  347.       ENDIF
  348.       ELSEIF (Nd.EQ.1) THEN
  349.       IF (I.NE.D(1)) THEN
  350.       Nu=Nu+1
  351.       U(Nu)=I
  352.       ENDIF
  353.       ENDIF
  354. 66    CONTINUE
  355.  
  356.       IF (Nu+Nd.NE.4) THEN
  357.       Print *,'Nd=,Nu=',Nd,Nu
  358.       ENDIF
  359.  
  360.       IF (Nu.EQ.2) THEN
  361.       IF ((U(1)+1).NE.U(2)) THEN
  362.       jj=U(1)
  363.       U(1)=U(2)
  364.       U(2)=jj
  365.       ENDIF
  366.       ENDIF
  367.  
  368.       N1 =U(1)
  369.       N2 =U(2)
  370.       N3 =D(1)
  371.       N4 =D(2)
  372.  
  373.       IF (Nd.EQ.2) THEN
  374.  
  375.       MINI2=0.D0
  376.       IF (0.D0.GE.KKK(N2))      MINI2=KKK(N2)
  377.  
  378.       MINI4=0.D0
  379.       IF (0.D0.GE.KKK(N4))      MINI4=KKK(N4)
  380.  
  381.       MAXI2=0.D0
  382.       IF (0.D0.LE.KKK(N2))      MAXI2=KKK(N2)
  383.  
  384.       MAXI4=0.D0
  385.       IF (0.D0.LE.KKK(N4))      MAXI4=KKK(N4)
  386.  
  387.       PhiN3 =( KKK(N1) + MINI2 + MINI4 )
  388.      *  *(TETAC(KP,N3)-TETAC(KP,N2))
  389.      *  +( MAXI4-MINI2 )*(TETAC(KP,N3)-TETAC(KP,N1))
  390.  
  391.       PhiN4 =( KKK(N1) + MINI2 + MINI4 )
  392.      * *(TETAC(KP,N4)-TETAC(KP,N1))
  393.      * +( MAXI2-MINI4 )*(TETAC(KP,N4)-TETAC(KP,N2))
  394.  
  395.       ZVGG(N3)=PhiN3
  396.       ZVGG(N4)=PhiN4
  397.  
  398.       r= -ZVGG(N3) / (ZVGG(N4)+XPETIT)
  399.       IF (r.GT.0.D0) THEN
  400.       IF (r.GT.1.D0) THEN
  401.       ZVGG(N3)=ZVGG(N3)+ZVGG(N4)
  402.       ZVGG(N4)=0.D0
  403.       ELSE
  404.       ZVGG(N4)=ZVGG(N3)+ZVGG(N4)
  405.       ZVGG(N3)=0.D0
  406.       ENDIF
  407.       ENDIF
  408.  
  409.       ELSEIF (ND.EQ.1) THEN
  410.  
  411.       PhiQ =0.5D0*(UMI(KP,1)*(n(1,1)+n(1,4)) +
  412.      * UMI(KP,2)*(n(2,1)+n(2,4)))*(TETAC(KP,3)-TETAC(KP,1))+
  413.      * 0.5D0*(UMI(KP,1)*(n(1,1)+n(1,2)) +
  414.      * UMI(KP,2)*(n(2,1)+n(2,2)))*(TETAC(KP,4)-TETAC(KP,2))
  415.  
  416.       ZVGG(N3) = PhiQ
  417.       ENDIF
  418.  
  419.       ELSEIF (NP.EQ.3) THEN
  420.  
  421. **************************************************
  422. *     Schéma PSI (Min-Mod) Pour des Triangles :
  423. **************************************************
  424.  
  425. *  U Noeuds amonts
  426. *  D Noeuds avals .
  427. **********************
  428.  
  429.       Nd=0
  430.       Nu=0
  431.  
  432.       DO 13 I=1,NP
  433.  
  434.       IF (KKK(I).GT.0.D0) THEN
  435.       Nd=Nd+1
  436.       D(Nd)=I
  437.       ELSE
  438.       Nu=Nu+1
  439.       U(Nu)=I
  440.       ENDIF
  441.  
  442. 13    CONTINUE
  443.  
  444.       IF (Nd.EQ.1) THEN
  445.  
  446. ***************************
  447. * 1 Target :
  448. ***************************
  449.  
  450.       N1 = D(1)
  451.       N2 = U(1)
  452.       N3 = U(2)
  453.  
  454.       ZVGG(N1) = PhiSource(KP) + KKK(1)*TETAC(KP,1)+
  455.      * KKK(2)*TETAC(KP,2) + KKK(3)*TETAC(KP,3)
  456.       ZVGG(N2) = 0.D0
  457.       ZVGG(N3) = 0.D0
  458.  
  459.       ELSEIF (Nd.EQ.2) THEN
  460.  
  461. ***************************
  462. * 2 Targets :
  463. ***************************
  464.  
  465.       N1 = D(1)
  466.       N2 = D(2)
  467.       N3 = U(1)
  468.  
  469.       ZVGG(N1) = KKK(N1)*(TETAC(KP,N1)-TETAC(KP,N3))-PhiSource(KP)*
  470.      * KKK(N1)/(KKK(N3)+XPETIT)
  471.  
  472.       ZVGG(N2) = KKK(N2)*(TETAC(KP,N2)-TETAC(KP,N3))-PhiSource(KP)*
  473.      * KKK(N2)/(KKK(N3)+XPETIT)
  474.  
  475.       ZVGG(N3) = 0.D0
  476.  
  477.       r= -ZVGG(N1) / (ZVGG(N2)+XPETIT)
  478.       IF (r.GT.0.D0) THEN
  479.       IF (r.GT.1.D0) THEN
  480.       ZVGG(N1)=ZVGG(N1)+ZVGG(N2)
  481.       ZVGG(N2)=0.D0
  482.       ELSE
  483.       ZVGG(N2)=ZVGG(N1)+ZVGG(N2)
  484.       ZVGG(N1)=0.D0
  485.       ENDIF
  486.       ENDIF
  487.  
  488.       ENDIF
  489.  
  490.       ENDIF
  491.  
  492.       DO 7 I=1,NP
  493.       SBF(KP,I)=ZVGG(I)
  494. 7     CONTINUE
  495.  
  496. 7008  CONTINUE
  497. **************************************************
  498.  
  499.       ELSEIF (IKOMP.EQ.1) THEN
  500.  
  501. * Prog à compléter
  502.       RETURN
  503.  
  504.       ENDIF
  505.  
  506. * Pas de temps DT :
  507.  
  508.       DT=100.
  509.  
  510.       DO 80066 K1=KDEB,KFIN
  511.       DO 32 J=1,NP
  512.       IP=IPADU(LE(J,K1))
  513.       DT1=1.D0 /( ZTE(IP)+XPETIT)
  514.  
  515. c     IF (Abs(DT1).GE.1000) THEN
  516. c     Print *,'DT1 trop grand =',DT1
  517. c     ENDIF
  518.  
  519.       IF (DT1.LE.DT) DT=DT1
  520. 32    CONTINUE
  521. 80066 CONTINUE
  522.  
  523.  
  524.       DTT1=DT
  525.       DTT2=DT
  526.       DTT3=0.D0
  527.       DIAEL=0.D0
  528.       NUEL=0
  529.  
  530. *      PRINT *, 'PAS DE TEMPS = ', DT
  531.  
  532.  
  533. ****************************************************
  534. * Contribution du terme diffusif:
  535.  
  536.       DO 80068 I=1,NP
  537.       DO 80067 J=1,NP
  538.       DO 8 K=KDEB,KFIN
  539.       KP=K-KDEB+1
  540.  
  541.       ZVGT=AIRE(KP)*(
  542.      &     HR(K,I,1)*HR(K,J,1)*ALFT(KP)
  543.      &+    HR(K,I,2)*HR(K,J,2)*ALFT(KP)
  544.      &+    (ALFT(KP)*AL2(KP)+ALFT(KP)*AH2(KP))/12.D0
  545.      &     *VGGT(J,I)*QUA4  )
  546.  
  547.  
  548.       SBF(KP,I)=SBF(KP,I)+ TETAD(KP,J)*ZVGT
  549.  
  550. 8     CONTINUE
  551. 80067 CONTINUE
  552. 80068 CONTINUE
  553. ****************************************************
  554.  
  555. C     Fin de l'accumulation dans SBF.
  556. C     On ajoute ces incréments G.
  557.  
  558.  
  559.       DO 80069 I=1,NP
  560.       DO 7017 K=KDEB,KFIN
  561.       KP=K-KDEB+1
  562.       NF=IPADL(LE(I,K))
  563.  
  564.       G(NF) = G(NF)+SBF(KP,I)
  565.  
  566. 7017  CONTINUE
  567. 80069 CONTINUE
  568.  
  569.  
  570. 7001  CONTINUE
  571.       IPAS=1
  572.  
  573.       RETURN
  574.  
  575.  
  576. C                           ********
  577. C                           *  3D  *
  578. C                           ********
  579.  
  580.  10   CONTINUE
  581.  
  582.       IF (IPAS.EQ.0) CALL CALHRH(QGGT,Q1,Q2,Q3,IES)
  583.       CUB8 = 0.D0
  584.       IF (NP.EQ.8) CUB8=1.D0
  585.  
  586. C
  587. C Calcul du nombre de paquets de LRV éléments
  588. C
  589.       NNN=MOD(NEL,LRV)
  590.       IF(NNN.EQ.0) NPACK=NEL/LRV
  591.       IF(NNN.NE.0) NPACK=1+(NEL-NNN)/LRV
  592.       KPACKD=1
  593.       KPACKF=NPACK
  594. C
  595. C ******* BOUCLE SUR LES PAQUETS DE LRV ELEMENTS **********
  596. C
  597.       DO 8001 KPACK=KPACKD,KPACKF
  598. C
  599. C ======= A L'INTERIEUR DE CHAQUE PAQUET DE LRV ELEMENTS =======
  600. C
  601. C 1. Calcul des limites du paquet courant.
  602.       KDEB=1+(KPACK-1)*LRV
  603.       KFIN=MIN(NEL,KDEB+LRV-1)
  604. C
  605.       DO 8002 K=KDEB,KFIN
  606.       KP=K-KDEB+1
  607.       NK=K+K0
  608.       NK1=(1-IND1)*(NK-1)+1
  609.       ALFT(KP)=ALFE(NK1)+XPETIT
  610.       AIRE(KP)=VOLU(NK)
  611.       AL(KP)=COTE(NK,1)
  612.       AH(KP)=COTE(NK,2)
  613.       AP(KP)=COTE(NK,3)
  614.  
  615.       CFM(KP)=AL(KP)*AH(KP)/AP(KP)+AL(KP)*AP(KP)/AH(KP)+
  616.      &        AP(KP)*AH(KP)/AL(KP)
  617. C     CF1(KP)=AL(KP)*AH(KP)/AP(KP)
  618. C     CF2(KP)=AL(KP)*AP(KP)/AH(KP)
  619. C     CF3(KP)=AP(KP)*AH(KP)/AL(KP)
  620.       XMH(KP)=(AL(KP)+AH(KP)+AP(KP))/3.D0
  621.       AL2(KP)=1.D0/AL(KP)/AL(KP)
  622.       AH2(KP)=1.D0/AH(KP)/AH(KP)
  623.       AP2(KP)=1.D0/AP(KP)/AP(KP)
  624.  
  625.  8002 CONTINUE
  626.       IF((IKOMP.EQ.0.AND.IKAS.EQ.5).OR.
  627.      &(IKOMP.EQ.1.AND.IKAS.EQ.6))THEN
  628.       DO 8003 K=KDEB,KFIN
  629.       KP=K-KDEB+1
  630.       NK=K+K0
  631.       ALFT(KP)=ALFT(KP)+ALT(NK)/SGT
  632.  8003 CONTINUE
  633.       ENDIF
  634.  
  635. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  636.  
  637. C Initialisation des UMI avant accumulation
  638.       DO 8005 K=KDEB,KFIN
  639.       KP=K-KDEB+1
  640.       UMI(KP,1)=XPETIT
  641.       UMI(KP,2)=XPETIT
  642.       UMI(KP,3)=XPETIT
  643.  8005 CONTINUE
  644.  
  645.       DO 80070 I=1,NP
  646.       DO 8006 K=KDEB,KFIN
  647.       KP=K-KDEB+1
  648.       NF=IPADL(LE(I,K))
  649.       NFU=IPADU(LE(I,K))
  650.       NFU=(1-INDU)*(NFU-1)+1
  651.       DR(KP,I)=DRR(I,K)
  652.       UMI(KP,1)=UMI(KP,1)+UN(NFU,1)*DR(KP,I)
  653.       UMI(KP,2)=UMI(KP,2)+UN(NFU,2)*DR(KP,I)
  654.       UMI(KP,3)=UMI(KP,3)+UN(NFU,3)*DR(KP,I)
  655.  
  656.       TETAC(KP,I)=HRN(NF)
  657.       TETAD(KP,I)=TN(NF)
  658. 8006  CONTINUE
  659. 80070 CONTINUE
  660.  
  661. C
  662. C Initialisation de la variable d'accumulation SBF au terme source
  663. C
  664. C       write(6,*)' IKomp,ikas=',IKomp,ikas
  665. C       write(6,*)' IKS,IND1,INDU=',IKS,IND1,INDU
  666.  
  667.       IF(IKOMP.EQ.0)THEN
  668.  
  669.          DO 80021 K=KDEB,KFIN
  670.          KP=K-KDEB+1
  671.          NK=K+K0
  672.          NKS=(1-IKS)*(NK-1)+1
  673.          QT(KP)=QE(NKS)
  674. 80021    CONTINUE
  675.  
  676.          DO 80071 I=1,NP
  677.          DO 80062 K=KDEB,KFIN
  678.          KP=K-KDEB+1
  679.          SBF(KP,I)=-QT(KP)*DR(KP,I)
  680. 80062    CONTINUE
  681. 80071    CONTINUE
  682.  
  683.       ELSEIF(IKOMP.EQ.1)THEN
  684.  
  685.          DO 80023 K=KDEB,KFIN
  686.          KP=K-KDEB+1
  687.          NK=K+K0
  688.          NKS=(1-IKS)*(NK-1)+1
  689.          QT(KP)=QE(NKS)
  690. 80023    CONTINUE
  691.  
  692.          DO 80072 I=1,NP
  693.          DO 80064 K=KDEB,KFIN
  694.          KP=K-KDEB+1
  695.          SBF(KP,I)=-QT(KP)*DR(KP,I)
  696. 80064    CONTINUE
  697. 80072    CONTINUE
  698.  
  699.       ENDIF
  700.  
  701.       DO 8007 K=KDEB,KFIN
  702.       KP=K-KDEB+1
  703.       UMI(KP,1)=UMI(KP,1)/AIRE(KP)
  704.       UMI(KP,2)=UMI(KP,2)/AIRE(KP)
  705.       UMI(KP,3)=UMI(KP,3)/AIRE(KP)
  706.       UM(KP)=UMI(KP,1)*UMI(KP,1)+UMI(KP,2)*UMI(KP,2)
  707.      &      +UMI(KP,3)*UMI(KP,3)
  708.       UM(KP)=SQRT(UM(KP))
  709. 8007  CONTINUE
  710.  
  711.  
  712. ********************************************************
  713. *          Debut de PSI:
  714.  
  715.  
  716.  
  717.  
  718.  
  719.  
  720.  
  721.  
  722. ********************************************************
  723.  
  724.  
  725. C Le coeur du calcul ...
  726.  
  727. *      IF(IKOMP.EQ.0)THEN
  728.  
  729. C      DO 8014 I=1,NP
  730. C      DO 8014 K=KDEB,KFIN
  731. C      KP=K-KDEB+1
  732.  
  733. *      ELSEIF(IKOMP.EQ.1)THEN
  734. *      RETURN
  735.  
  736. *      ENDIF
  737.  
  738. C
  739. C     Fin de l'accumulation dans SBF.
  740. C     On ajoute ces incréments G.
  741. C
  742.       DO 80073 I=1,NP
  743.       DO 8017 K=KDEB,KFIN
  744.  
  745.       KP=K-KDEB+1
  746.       NF=IPADL(LE(I,K))
  747.       G(NF) = G(NF)+SBF(KP,I)
  748. 8017  CONTINUE
  749. 80073 CONTINUE
  750.  
  751. 8001  CONTINUE
  752.  
  753.  
  754.       IPAS=1
  755.       RETURN
  756. 1002  FORMAT(10(1X,1PE11.4))
  757.  
  758.       END
  759.  
  760.  
  761.  
  762.  
  763.  
  764.  
  765.  

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