Télécharger la plaquette Cast3M


fr en

Télécharger yctscl.eso

Retour à la liste

Numérotation des lignes : 

yctscl

  1. C YCTSCL    SOURCE    CHAT      05/01/13    04:16:35     5004
  2.       SUBROUTINE YCTSCL
  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,IES,NP,IAXI,
  12.      & IPADI,IPADS,IPADF,IKOMP,IKAS,
  13.      & ALFE,IND1,UN,INDU,NPTS,TN,NPTD,QE,IKS,
  14.      & HRN,G,NPTI,
  15.      & ALT,SGT,
  16.      & VOLU,COTE,NELZ,IDCEN,IPG,
  17.      & DTM1,DT,DTT1,DTT2,NUEL,DIAEL,FN)
  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. C
  39. C Longueur des registres vectoriels de la machine cible
  40. C On prend 64 pour ne pas augmenter la taille des tableaux
  41. C nécessaires à la vectorisation.
  42. C
  43.       PARAMETER(LRV=64)
  44.  
  45.       DIMENSION UN(NPTS,IES),HRN(NPTI),TN(NPTD)
  46.       DIMENSION COTE(NELZ,IES),VOLU(NELZ),QE(*)
  47.       DIMENSION ALFE(*),ALT(*)
  48.  
  49.       DIMENSION IPADI(*),IPADS(*),IPADF(*),LE(NP,1)
  50.       DIMENSION HR(NEL,NP,IES),RPG(1),DRR(NP,NEL)
  51.  
  52.       DIMENSION BF(9,9)
  53.       DIMENSION QGGT(8,8),Q1(8,8),Q2(8,8),Q3(8,8)
  54.  
  55.       DIMENSION AIRE(LRV)
  56.       DIMENSION AL(LRV),AH(LRV),AP(LRV)
  57.       DIMENSION ALFT(LRV),QT(LRV)
  58.       DIMENSION UIX(LRV,9),UIY(LRV,9),UIZ(LRV,9)
  59.       DIMENSION TETAC(LRV,9),TETAD(LRV,9),TETA(LRV,9)
  60.       DIMENSION UMI(LRV,3)
  61.       DIMENSION SBF(LRV,9)
  62.       DIMENSION WT(LRV,9),CHGLD(LRV),CHGLP(LRV)
  63.  
  64.       REAL*8    G(NPTS),FN(NP,*)
  65.  
  66.       SAVE IPAS,QGGT,Q1,Q2,Q3
  67.  
  68.       DATA CD/1.D0/
  69.  
  70.       DATA IPAS/0/
  71. C************************************************************************
  72. C
  73. C           INITIALISATIONS DIVERSES
  74. C
  75.        ZERMA=XPETIT
  76.  
  77.       NK=K0
  78. C                           ********
  79. C                           *  2D  *
  80. C                           ********
  81.  
  82.  
  83.       IF(IES.EQ.3)GO TO 10
  84.  
  85.       IAX1=0
  86.       IAX2=0
  87.       IF(IAXI.EQ.1)IAX2=1
  88.       IF(IAXI.EQ.2)IAX1=1
  89.  
  90. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  91.  
  92. C                           DIFFERENCES TRIANGLE / QUADRANGLE
  93.       QUA4=0.D0
  94.       IF(NP.EQ.4)QUA4=1.D0
  95. C
  96. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  97.  
  98. C
  99. C Calcul du nombre de paquets de LRV éléments
  100. C
  101.       NNN=MOD(NEL,LRV)
  102.       IF(NNN.EQ.0) NPACK=NEL/LRV
  103.       IF(NNN.NE.0) NPACK=1+(NEL-NNN)/LRV
  104.       KPACKD=1
  105.       KPACKF=NPACK
  106. C
  107. C ******* BOUCLE SUR LES PAQUETS DE LRV ELEMENTS **********
  108. C
  109.       DO 7001 KPACK=KPACKD,KPACKF
  110. C
  111. C ======= A L'INTERIEUR DE CHAQUE PAQUET DE LRV ELEMENTS =======
  112. C
  113. C 1. Calcul des limites du paquet courant.
  114.       KDEB=1+(KPACK-1)*LRV
  115.       KFIN=MIN(NEL,KDEB+LRV-1)
  116. C
  117.       DO 7002 K=KDEB,KFIN
  118.       KP=K-KDEB+1
  119.       NK=K+K0
  120.       NK1=(1-IND1)*(NK-1)+1
  121.       ALFT(KP)=ALFE(NK1)+ZERMA
  122.       AIRE(KP)=VOLU(NK)
  123.       AL(KP)=COTE(NK,1)
  124.       AH(KP)=COTE(NK,2)
  125.  7002 CONTINUE
  126.  
  127.       IF((IKOMP.EQ.0.AND.IKAS.EQ.5).OR.
  128.      &(IKOMP.EQ.1.AND.IKAS.EQ.6))THEN
  129.       DO 7003 K=KDEB,KFIN
  130.       KP=K-KDEB+1
  131.       NK=K+K0
  132.       ALFT(KP)=ALFT(KP)+ALT(NK)/SGT
  133.  7003 CONTINUE
  134.       ENDIF
  135.  
  136. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  137.  
  138.       DO 7006 I=1,NP
  139.       DO 7006 K=KDEB,KFIN
  140.       KP=K-KDEB+1
  141.       NF=IPADF(LE(I,K))
  142.       NU=IPADS(LE(I,K))
  143.       NI=IPADI(LE(I,K))
  144.       NFU=(1-INDU)*(NU-1)+1
  145.       UIX(KP,I)=UN(NFU,1)
  146.       UIY(KP,I)=UN(NFU,2)
  147.       TETAC(KP,I)=HRN(NI)
  148.       TETAD(KP,I)=TN(NF)
  149.  7006 CONTINUE
  150.  
  151.       CALL KSUPG1(WT,UMI,CHGLD,CHGLP,KDEB,KFIN,LRV,
  152.      &HRN,IPADI,UN,ALFT,NPTS,NEL,NP,DRR,HR,FN,
  153.      &AIRE,AL,AH,AP,IDCEN,IPADS,LE,QUA4,IKOMP,
  154.      &DTM1,DT,DTT1,DTT2,DIAEL,NUEL)
  155. C
  156. C Initialisation de la variable d'accumulation SBF au terme source
  157. C
  158.  
  159.       IF(IKOMP.EQ.0)THEN
  160.  
  161.          DO 70021 K=KDEB,KFIN
  162.          KP=K-KDEB+1
  163.          NK=K+K0
  164.          NKS=(1-IKS)*(NK-1)+1
  165.          QT(KP)=QE(NKS)
  166. 70021    CONTINUE
  167.  
  168.          IF(IPG.EQ.0)THEN
  169.          DO 70062 I=1,NP
  170.          DO 70062 K=KDEB,KFIN
  171.          KP=K-KDEB+1
  172.          SBF(KP,I)=-QT(KP)*DRR(I,K)
  173. 70062    CONTINUE
  174.          ELSE
  175.          DO 71062 I=1,NP
  176.          DO 71062 K=KDEB,KFIN
  177.          KP=K-KDEB+1
  178.          SBF(KP,I)=-QT(KP)*WT(KP,I)
  179. 71062    CONTINUE
  180.          ENDIF
  181.  
  182.       ELSEIF(IKOMP.EQ.1)THEN
  183.  
  184.          DO 70023 K=KDEB,KFIN
  185.          KP=K-KDEB+1
  186.          NK=K+K0
  187.          NKS=(1-IKS)*(NK-1)+1
  188.          QT(KP)=QE(NKS)
  189. 70023    CONTINUE
  190.  
  191.        IF(IPG.EQ.0)THEN
  192.          DO 70064 I=1,NP
  193.          DO 70064 K=KDEB,KFIN
  194.          KP=K-KDEB+1
  195.          SBF(KP,I)=-QT(KP)*DRR(I,K)
  196. 70064    CONTINUE
  197.        ELSE
  198.          DO 71064 I=1,NP
  199.          DO 71064 K=KDEB,KFIN
  200.          KP=K-KDEB+1
  201.          SBF(KP,I)=-QT(KP)*WT(KP,I)
  202. 71064    CONTINUE
  203.        ENDIF
  204.  
  205.       ENDIF
  206.  
  207. C Le coeur du calcul ...
  208.  
  209.       IF(IKOMP.EQ.0)THEN
  210.  
  211.       DO 7014 I=1,NP
  212.       DO 7014 J= 1,NP
  213.       DO 7014 K=KDEB,KFIN
  214.       KP=K-KDEB+1
  215.       ZVGG=AIRE(KP)*CHGLP(KP)*VGGT(J,I)
  216.  
  217.       ZVGT=AIRE(KP)*(
  218.      &     HR(K,I,1)*HR(K,J,1)*ALFT(KP)
  219.      &+    HR(K,I,2)*HR(K,J,2)*ALFT(KP)
  220.      &+    CHGLD(KP)*VGGT(J,I) )
  221.  
  222.       V2=(UMI(KP,1)*HR(K,J,1)+UMI(KP,2)*HR(K,J,2))*WT(KP,I)
  223.  
  224.       SBF(KP,I)=SBF(KP,I)+TETAC(KP,J)*(ZVGG+V2)+ TETAD(KP,J)*ZVGT
  225.  
  226.  7014 CONTINUE
  227.  
  228.       ELSEIF(IKOMP.EQ.1)THEN
  229.  
  230.       DO 7015 I=1,NP
  231.       DO 7015 J= 1,NP
  232.       DO 7015 K=KDEB,KFIN
  233.       KP=K-KDEB+1
  234.       ZVGG=AIRE(KP)*CHGLP(KP)*VGGT(J,I)
  235.  
  236.       ZVGT=AIRE(KP)*(
  237.      &     HR(K,I,1)*HR(K,J,1)*ALFT(KP)
  238.      &+    HR(K,I,2)*HR(K,J,2)*ALFT(KP)
  239.      &+    CHGLD(KP)*VGGT(J,I) )
  240.  
  241.       V2=(UIX(KP,J)*HR(K,J,1)+UIY(KP,J)*HR(K,J,2))*WT(KP,I)
  242.  
  243.       SBF(KP,I)=SBF(KP,I)+TETAC(KP,J)*(ZVGG+V2)+ TETAD(KP,J)*ZVGT
  244.  
  245.  7015 CONTINUE
  246.  
  247.       ENDIF
  248. C
  249. C     Fin de l'accumulation dans SBF.
  250. C     On ajoute ces incréments G.
  251. C
  252.       DO 7017 I=1,NP
  253.       DO 7017 K=KDEB,KFIN
  254.       KP=K-KDEB+1
  255.       NF=IPADS(LE(I,K))
  256.       G(NF) = G(NF)-SBF(KP,I)
  257.  7017 CONTINUE
  258.  
  259.  7001 CONTINUE
  260.  
  261. C     WRITE(6,*)' G DANS YCTSCL '
  262. C     WRITE(6,1984)(M,G(M),M=1,NPTS)
  263.  1984 FORMAT(7(1X,I4,2X,1PE11.4))
  264.  
  265. C     CALL ARRET(0)
  266.       IPAS=1
  267.       RETURN
  268.  
  269. C                           ********
  270. C                           *  3D  *
  271. C                           ********
  272.  
  273.  10   CONTINUE
  274.  
  275.       IF(IPAS.EQ.0)CALL CALHRH(QGGT,Q1,Q2,Q3,IES)
  276.       CUB8=0.D0
  277.       IF(NP.EQ.8)CUB8=1.D0
  278.  
  279. C
  280. C Calcul du nombre de paquets de LRV éléments
  281. C
  282.       NNN=MOD(NEL,LRV)
  283.       IF(NNN.EQ.0) NPACK=NEL/LRV
  284.       IF(NNN.NE.0) NPACK=1+(NEL-NNN)/LRV
  285.       KPACKD=1
  286.       KPACKF=NPACK
  287. C
  288. C ******* BOUCLE SUR LES PAQUETS DE LRV ELEMENTS **********
  289. C
  290.       DO 8001 KPACK=KPACKD,KPACKF
  291. C
  292. C ======= A L'INTERIEUR DE CHAQUE PAQUET DE LRV ELEMENTS =======
  293. C
  294. C 1. Calcul des limites du paquet courant.
  295.       KDEB=1+(KPACK-1)*LRV
  296.       KFIN=MIN(NEL,KDEB+LRV-1)
  297. C
  298.       DO 8002 K=KDEB,KFIN
  299.       KP=K-KDEB+1
  300.       NK=K+K0
  301.       NK1=(1-IND1)*(NK-1)+1
  302.       ALFT(KP)=ALFE(NK1)+ZERMA
  303.       AIRE(KP)=VOLU(NK)
  304.       AL(KP)=COTE(NK,1)
  305.       AH(KP)=COTE(NK,2)
  306.       AP(KP)=COTE(NK,3)
  307.  8002 CONTINUE
  308.       IF((IKOMP.EQ.0.AND.IKAS.EQ.5).OR.
  309.      &(IKOMP.EQ.1.AND.IKAS.EQ.6))THEN
  310.       DO 8003 K=KDEB,KFIN
  311.       KP=K-KDEB+1
  312.       NK=K+K0
  313.       ALFT(KP)=ALFT(KP)+ALT(NK)/SGT
  314.  8003 CONTINUE
  315.       ENDIF
  316.  
  317. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  318.  
  319.       DO 8006 I=1,NP
  320.       DO 8006 K=KDEB,KFIN
  321.       KP=K-KDEB+1
  322.       NF=IPADF(LE(I,K))
  323.       NU=IPADS(LE(I,K))
  324.       NI=IPADI(LE(I,K))
  325.       NFU=(1-INDU)*(NU-1)+1
  326.       UIX(KP,I)=UN(NFU,1)
  327.       UIY(KP,I)=UN(NFU,2)
  328.       UIZ(KP,I)=UN(NFU,3)
  329.       TETAC(KP,I)=HRN(NI)
  330.       TETAD(KP,I)=TN(NF)
  331.  8006 CONTINUE
  332.  
  333.       CALL KSUPG1(WT,UMI,CHGLD,CHGLP,KDEB,KFIN,LRV,
  334.      &HRN,IPADI,UN,ALFT,NPTS,NEL,NP,DRR,HR,FN,
  335.      &AIRE,AL,AH,AP,IDCEN,IPADS,LE,CUB8,IKOMP,
  336.      &DTM1,DT,DTT1,DTT2,DIAEL,NUEL)
  337.  
  338. C
  339. C Initialisation de la variable d'accumulation SBF au terme source
  340. C                                                                      M
  341.       IF(IKOMP.EQ.0)THEN
  342.  
  343.          DO 80021 K=KDEB,KFIN
  344.          KP=K-KDEB+1
  345.          NK=K+K0
  346.          NKS=(1-IKS)*(NK-1)+1
  347.          QT(KP)=QE(NKS)
  348. 80021    CONTINUE
  349.  
  350.        IF(IPG.EQ.0)THEN
  351.          DO 80062 I=1,NP
  352.          DO 80062 K=KDEB,KFIN
  353.          KP=K-KDEB+1
  354.          SBF(KP,I)=-QT(KP)*DRR(I,K)
  355. 80062    CONTINUE
  356.        ELSE
  357.          DO 81062 I=1,NP
  358.          DO 81062 K=KDEB,KFIN
  359.          KP=K-KDEB+1
  360.          SBF(KP,I)=-QT(KP)*WT(KP,I)
  361. 81062    CONTINUE
  362.        ENDIF
  363.  
  364.       ELSEIF(IKOMP.EQ.1)THEN
  365.  
  366.          DO 80023 K=KDEB,KFIN
  367.          KP=K-KDEB+1
  368.          NK=K+K0
  369.          NKS=(1-IKS)*(NK-1)+1
  370.          QT(KP)=QE(NKS)
  371. 80023    CONTINUE
  372.  
  373.        IF(IPG.EQ.0)THEN
  374.          DO 80064 I=1,NP
  375.          DO 80064 K=KDEB,KFIN
  376.          KP=K-KDEB+1
  377.          SBF(KP,I)=-QT(KP)*DRR(I,K)
  378. 80064    CONTINUE
  379.        ELSE
  380.          DO 81064 I=1,NP
  381.          DO 81064 K=KDEB,KFIN
  382.          KP=K-KDEB+1
  383.          SBF(KP,I)=-QT(KP)*WT(KP,I)
  384. 81064    CONTINUE
  385.        ENDIF
  386.  
  387.       ENDIF
  388.  
  389. C Le coeur du calcul ...
  390.  
  391.       IF(IKOMP.EQ.0)THEN
  392.  
  393.       DO 8014 I=1,NP
  394.       DO 8014 J= 1,NP
  395.       DO 8014 K=KDEB,KFIN
  396.       KP=K-KDEB+1
  397.  
  398.       ZVGG=AIRE(KP)*CHGLP(KP)*QGGT(J,I)
  399.  
  400.       ZVGT=AIRE(KP)*(
  401.      &     HR(K,I,1)*HR(K,J,1)*ALFT(KP)
  402.      &+    HR(K,I,2)*HR(K,J,2)*ALFT(KP)
  403.      &+    HR(K,I,3)*HR(K,J,3)*ALFT(KP) )
  404.      &+    CHGLD(KP)*QGGT(J,I)
  405.  
  406.       V2=UMI(KP,1)*HR(K,J,1)+UMI(KP,2)*HR(K,J,2)+UMI(KP,3)*HR(K,J,3)
  407.  
  408.       SBF(KP,I)=SBF(KP,I)
  409.      &         +TETAC(KP,J)*(ZVGG+V2*WT(KP,I))+TETAD(KP,J)*ZVGT
  410.  
  411.  8014 CONTINUE
  412.  
  413.       ELSEIF(IKOMP.EQ.1)THEN
  414.  
  415.       DO 8015 I=1,NP
  416.       DO 8015 J= 1,NP
  417.       DO 8015 K=KDEB,KFIN
  418.  
  419.       KP=K-KDEB+1
  420.  
  421.       ZVGG=AIRE(KP)*CHGLP(KP)*QGGT(J,I)
  422.  
  423.       ZVGT=AIRE(KP)*(
  424.      &     HR(K,I,1)*HR(K,J,1)*ALFT(KP)
  425.      &+    HR(K,I,2)*HR(K,J,2)*ALFT(KP)
  426.      &+    HR(K,I,3)*HR(K,J,3)*ALFT(KP) )
  427.      &+    CHGLD(KP)*QGGT(J,I)
  428.  
  429.       V2=UIX(KP,J)*HR(K,J,1)+UIY(KP,J)*HR(K,J,2)+UIZ(KP,J)*HR(K,J,3)
  430.  
  431.       SBF(KP,I)=SBF(KP,I)
  432.      &         +TETAC(KP,J)*(ZVGG+V2*WT(KP,I))+TETAD(KP,J)*ZVGT
  433.  
  434.  8015 CONTINUE
  435.  
  436.       ENDIF
  437.  
  438. C
  439. C     Fin de l'accumulation dans SBF.
  440. C     On ajoute ces incréments G.
  441. C
  442.       DO 8017 I=1,NP
  443.       DO 8017 K=KDEB,KFIN
  444.       KP=K-KDEB+1
  445.       NF=IPADS(LE(I,K))
  446.       G(NF) = G(NF)-SBF(KP,I)
  447.  8017 CONTINUE
  448.  
  449.  8001 CONTINUE
  450.  
  451. C     WRITE(6,*)' G DANS YCTSCL '
  452. C     WRITE(6,1984)(M,G(M),M=1,NPTS)
  453.  
  454. C     CALL ARRET(0)
  455.       IPAS=1
  456.       RETURN
  457.  1002 FORMAT(10(1X,1PE11.4))
  458.       END
  459.  
  460.  
  461.  
  462.  
  463.  
  464.  
  465.  
  466.  
  467.  

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