Télécharger ksupg1.eso

Retour à la liste

Numérotation des lignes :

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

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