Télécharger la plaquette Cast3M


fr en

Télécharger zctscl.eso

Retour à la liste

Numérotation des lignes : 

zctscl

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

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