Télécharger cli281.eso

Retour à la liste

Numérotation des lignes :

  1. C CLI281 SOURCE KK2000 14/04/10 21:15:04 8032
  2. SUBROUTINE CLI281(NSP,MELEMF,MELEMC,MELECB,MELEFC,INORM,ICHPVO,
  3. & ICHPSU,LRECP,LRECV,IROC,IVITC,IPC,IYN,ICHLIM,ICHRES,ICHRLI)
  4. C************************************************************************
  5. C
  6. C PROJET : CASTEM 2000
  7. C
  8. C NOM : CLI281
  9. C
  10. C DESCRIPTION : Subroutine appellée par CLIM22
  11. C calcul de RESIDU et CLIM at the board
  12. C OPTION: 'INOU' 2D
  13. C
  14. C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
  15. C
  16. C AUTEUR : S.Kudriakov, DEN/DM2S/SFME/LTMF
  17. C
  18. C************************************************************************
  19. C
  20. C APPELES (Calcul) :
  21. C
  22. C************************************************************************
  23. C
  24. C HISTORIQUE (Anomalies et modifications éventuelles)
  25. C
  26. C HISTORIQUE :
  27. C
  28. C************************************************************************
  29. C
  30. IMPLICIT INTEGER(I-N)
  31. IMPLICIT REAL*8(A-H,O-Z)
  32.  
  33. -INC CCOPTIO
  34. -INC SMLMOTS
  35. -INC SMELEME
  36. POINTEUR MELEFC.MELEME
  37. -INC SMLENTI
  38. POINTEUR MLEMC.MLENTI, MLEMCB.MLENTI,MLEMF.MLENTI
  39. -INC SMCHPOI
  40. POINTEUR MPNORM.MPOVAL, MPVOL.MPOVAL, MPSURF.MPOVAL, MPRC.MPOVAL,
  41. & MPVC.MPOVAL, MPPC.MPOVAL, MPYN.MPOVAL, MPLIM.MPOVAL,
  42. & MPRES.MPOVAL, MPRLI.MPOVAL
  43.  
  44. INTEGER MELEMF,MELEMC,MELECB,INORM,ICHPVO,ICHPSU, IROC,IVITC,IPC
  45. & ,IYN,ICHLIM,ICHRES,ICHRLI,ICEL,NFAC,IFAC
  46. & ,NGF,NGC,NLF,NLC,NLCB,LRECP,LRECV,I,NSP,NESP
  47. REAL*8 VOLU,SURF,RC,PC,UXC,UYC,UZC,GAMC,CNX,CNY,CNZ,CTX,CTY,CTZ
  48. & ,CT2X,CT2Y,CT2Z,RF,PF,UXF,UYF,UZF,TOP,BOT
  49. & ,UNC,UNF,UTF,UT2F,SF,GAMF,ECIN,PSRF,HTF,GM1
  50. & ,CELLT,UT2C,UTC,TEMP
  51. CHARACTER*(8) TYPE
  52. C------------------------------------------------------------
  53. -INC SMLREEL
  54. POINTEUR MLRECP.MLREEL, MLRECV.MLREEL
  55. C-------------------------------------------------------
  56. C********** Les CP's and CV's ***********************
  57. C-------------------------------------------------------
  58. SEGMENT GCONST
  59. REAL*8 GC(NSP)
  60. ENDSEGMENT
  61. POINTEUR CP.GCONST, CV.GCONST
  62. C-------------------------------------------------------------
  63. C******* Les fractionines massiques **************************
  64. C-------------------------------------------------------------
  65. SEGMENT FRAMAS
  66. REAL*8 YET(NSP)
  67. ENDSEGMENT
  68. POINTEUR YC.FRAMAS, YF.FRAMAS
  69. C-------------------------------------------------------------
  70. C********** Segments for the flux-vector *******************
  71. C-------------------------------------------------------------
  72. SEGMENT FUNEL
  73. REAL*8 FU(4+NSP)
  74. ENDSEGMENT
  75. POINTEUR flux2D.funel, flux3D.funel
  76. SEGINI FLUX2D
  77. SEGINI FLUX3D
  78. C------------------------------------------------------
  79. C**** KRIPAD pour la correspondance global/local
  80. C------------------------------------------------------
  81. CALL KRIPAD(MELEMC,MLEMC)
  82. CALL KRIPAD(MELECB,MLEMCB)
  83. CALL KRIPAD(MELEMF,MLEMF)
  84. C------------------------------------------------------
  85. C**** CHPOINTs de la table DOMAINE
  86. C------------------------------------------------------
  87. CALL LICHT(INORM,MPNORM,TYPE,ICEL)
  88. CALL LICHT(ICHPVO,MPVOL,TYPE,ICEL)
  89. CALL LICHT(ICHPSU,MPSURF,TYPE,ICEL)
  90. C------------------------------------------------------
  91. C**** CHPOINTs des variables
  92. C------------------------------------------------------
  93. CALL LICHT(IROC,MPRC,TYPE,ICEL)
  94. CALL LICHT(IVITC,MPVC,TYPE,ICEL)
  95. CALL LICHT(IPC,MPPC,TYPE,ICEL)
  96. CALL LICHT(IYN,MPYN,TYPE,ICEL)
  97. CALL LICHT(ICHLIM,MPLIM,TYPE,ICEL)
  98. CALL LICHT(ICHRES,MPRES,TYPE,ICEL)
  99. CALL LICHT(ICHRLI,MPRLI,TYPE,ICEL)
  100. C---------------------------------------------------------
  101. C**** Boucle sur le face pour le calcul des invariants de
  102. C Riemann et du flux
  103. C---------------------------------------------------------
  104. SEGACT MELEFC
  105. NFAC=MELEFC.NUM(/2)
  106. c write(*,*) 'hell', NFAC
  107. UZC=0.0D0
  108. CNZ=0.0D0
  109. CTZ=0.0D0
  110. CT2X=0.0D0
  111. CT2Y=0.0D0
  112. CT2Z=0.0D0
  113. DO 1 IFAC=1,NFAC,1
  114. NGF=MELEFC.NUM(1,IFAC)
  115. NGC=MELEFC.NUM(2,IFAC)
  116. NLF=MLEMF.LECT(NGF)
  117. NLC=MLEMC.LECT(NGC)
  118. NLCB=MLEMCB.LECT(NGF)
  119. VOLU=MPVOL.VPOCHA(NLC,1)
  120. SURF=MPSURF.VPOCHA(NLF,1)
  121. C----------------------------------------------
  122. C In CASTEM les normales sont sortantes
  123. C----------------------------------------------
  124. CNX=MPNORM.VPOCHA(NLF,1)
  125. CNY=MPNORM.VPOCHA(NLF,2)
  126. IF(IDIM.EQ.2)THEN
  127. CTX=-1.0D0*CNY
  128. CTY=CNX
  129. ELSE
  130. CNZ=MPNORM.VPOCHA(NLF,3)
  131. CTX=MPNORM.VPOCHA(NLF,4)
  132. CTY=MPNORM.VPOCHA(NLF,5)
  133. CTZ=MPNORM.VPOCHA(NLF,6)
  134. CT2X=MPNORM.VPOCHA(NLF,7)
  135. CT2Y=MPNORM.VPOCHA(NLF,8)
  136. CT2Z=MPNORM.VPOCHA(NLF,9)
  137. ENDIF
  138. C----------------------------------------
  139. SEGINI CP, CV
  140. MLRECP = LRECP
  141. MLRECV = LRECV
  142. SEGACT MLRECP, MLRECV
  143. DO 10 I=1,(NSP-1)
  144. CP.GC(I)=MLRECP.PROG(I)
  145. CV.GC(I)=MLRECV.PROG(I)
  146. 10 CONTINUE
  147. CP.GC(NSP)=MLRECP.PROG(NSP)
  148. CV.GC(NSP)=MLRECV.PROG(NSP)
  149. C----------------------------
  150. C Variables au centre
  151. C----------------------------
  152. RC=MPRC.VPOCHA(NLC,1)
  153. PC=MPPC.VPOCHA(NLC,1)
  154. UXC=MPVC.VPOCHA(NLC,1)
  155. UYC=MPVC.VPOCHA(NLC,2)
  156. IF(IDIM.EQ.3)UZC=MPVC.VPOCHA(NLC,3)
  157. SEGINI YC
  158. SEGACT MPYN
  159. DO 100 I=1,(NSP-1)
  160. YC.YET(I)=MPYN.VPOCHA(NLC,I)
  161. 100 CONTINUE
  162. c-------------------------------------------------------------
  163. c Computing GAMMA at the cell-center
  164. c-------------------------------------------------------------
  165. top=0.0D0
  166. bot=0.0D0
  167. do 102 i=1,(nsp-1)
  168. top=top+yc.yet(i)*(cp.gc(i)-cp.gc(nsp))
  169. bot=bot+yc.yet(i)*(cv.gc(i)-cv.gc(nsp))
  170. 102 continue
  171. top=cp.gc(nsp)+top
  172. bot=cv.gc(nsp)+bot
  173. GAMC=top/bot
  174. C-------------------------------------------------------------
  175. C We decide whether the boundary conditions are 'inlet' or
  176. C 'outlet' by looking at the UNC (UNC > 0.0 -- 'outlet')
  177. C-------------------------------------------------------------
  178. UNC=(UXC*CNX)+(UYC*CNY)+(UZC*CNZ)
  179. c write(*,*) 'before', UNC, NLC,NGC
  180. IF(UNC .GE. 0.0D0) THEN
  181. C--------------------------------------
  182. C 'OUTP' -- outlet with given pressure
  183. C--------------------------------------
  184. C-----------------------------------------
  185. C Variables à la face
  186. C-----------------------------------------
  187. PF=MPLIM.VPOCHA(NLCB,NSP+2)
  188. C---------------------------------------
  189. C******* On calcule UN, UT, UT2
  190. C---------------------------------------
  191. UTC=(UXC*CTX)+(UYC*CTY)+(UZC*CTZ)
  192. UT2C=(UXC*CT2X)+(UYC*CT2Y)+(UZC*CT2Z)
  193. C-----------------------------------------------
  194. C******* Densite, vitesse, pression sur le bord
  195. C-----------------------------------------------
  196. MPRLI.VPOCHA(NLCB,1)=RC
  197. MPRLI.VPOCHA(NLCB,2)=UXC
  198. MPRLI.VPOCHA(NLCB,3)=UYC
  199. IF(IDIM.EQ.3) MPRLI.VPOCHA(NLCB,4)=UZC
  200. MPRLI.VPOCHA(NLCB,IDIM+2)=PF
  201. do 104 i=1,(nsp-1)
  202. MPRLI.VPOCHA(NLCB,IDIM+2+I)=YC.YET(I)
  203. 104 continue
  204. C---------------------------------------------------
  205. C******* Probleme de Riemann entre l'etat gauche
  206. C RC,UNC,UTC,UT2C,PC et l'etat droite
  207. C RC,UNC,UTC,UT2C,PF
  208. C On utilise AUSM+
  209. C Flux dans le repaire normale
  210. C---------------------------------------------------
  211. NESP=NSP-1
  212. IF(IDIM.EQ.2)THEN
  213. CALL FAUSMP(NESP,
  214. & GAMC,RC,PC,UNC,UTC,
  215. & GAMC,RC,PF,UNC,UTC,
  216. & YC.YET,YC.YET,
  217. & FLUX2D.FU,
  218. & CELLT)
  219. C-------------------------------------------------------
  220. C******* Residuum (son SPG a le meme ordre que MELEFC)
  221. C-------------------------------------------------------
  222. MPRES.VPOCHA(IFAC,1)=-1*FLUX2D.FU(1)*SURF/VOLU
  223. MPRES.VPOCHA(IFAC,2)=-1*((FLUX2D.FU(2)*CNX)+
  224. & (FLUX2D.FU(3)*CTX))*SURF/VOLU
  225. MPRES.VPOCHA(IFAC,3)=-1*((FLUX2D.FU(2)*CNY)+
  226. & (FLUX2D.FU(3)*CTY))*SURF/VOLU
  227. MPRES.VPOCHA(IFAC,4)=-1*FLUX2D.FU(4)*SURF/VOLU
  228. do 105 i=1,(nsp-1)
  229. MPRES.VPOCHA(IFAC,4+I)=-1*FLUX2D.FU(4+I)*SURF/VOLU
  230. 105 continue
  231. ELSE
  232. CALL FAUSM3(NESP,
  233. & GAMC,RC,PC,UNC,UTC,UT2C,
  234. & GAMC,RC,PF,UNC,UTC,UT2C,
  235. & YC.YET,YC.YET,
  236. & FLUX3D.FU,
  237. & CELLT)
  238. C------------------------------------------------------
  239. C******* Residuum (son SPG a le meme ordre que MELEFC)
  240. C------------------------------------------------------
  241. MPRES.VPOCHA(IFAC,1)=-1*FLUX3D.FU(1)*SURF/VOLU
  242. MPRES.VPOCHA(IFAC,2)=-1*((FLUX3D.FU(2)*CNX)+
  243. & (FLUX3D.FU(3)*CTX)+(FLUX3D.FU(4)*CT2X))*SURF/VOLU
  244. MPRES.VPOCHA(IFAC,3)=-1*((FLUX3D.FU(2)*CNY)+
  245. & (FLUX3D.FU(3)*CTY)+(FLUX3D.FU(4)*CT2Z))*SURF/VOLU
  246. MPRES.VPOCHA(IFAC,4)=-1*((FLUX3D.FU(2)*CNZ)+
  247. & (FLUX3D.FU(3)*CTZ)+(FLUX3D.FU(4)*CT2Z))*SURF/VOLU
  248. MPRES.VPOCHA(IFAC,5)=-1*FLUX3D.FU(5)*SURF/VOLU
  249. do 106 i=1,(nsp-1)
  250. MPRES.VPOCHA(IFAC,5+I)=-1*FLUX3D.FU(5+I)*SURF/VOLU
  251. 106 continue
  252. ENDIF
  253. TEMP=-1000.0
  254. ELSE
  255. c write(*,*) 'insuuuuuuuu', NLCB
  256. C-------------------------------------------------------
  257. C 'INSU' -- 'inlet' b.c. with given Ht and S (entropy)
  258. C-------------------------------------------------------
  259. UZF=0.0D0
  260. UT2F=0.0D0
  261. C----------------------------
  262. C Variables à la face
  263. C----------------------------
  264. HTF=MPLIM.VPOCHA(NLCB,1)
  265. SF=MPLIM.VPOCHA(NLCB,2)
  266. SEGINI YF
  267. TEMP=1000.0
  268. DO 101 I=1,(NSP-1)
  269. YF.YET(I)=MPLIM.VPOCHA(NLCB,2+I)
  270. 101 CONTINUE
  271. UTF=0.0D0
  272. c-------------------------------------------------------------
  273. c Computing GAMMA at the face-center
  274. c-------------------------------------------------------------
  275. top=0.0D0
  276. bot=0.0D0
  277. do 103 i=1,(nsp-1)
  278. top=top+yf.yet(i)*(cp.gc(i)-cp.gc(nsp))
  279. bot=bot+yf.yet(i)*(cv.gc(i)-cv.gc(nsp))
  280. 103 continue
  281. top=cp.gc(nsp)+top
  282. bot=cv.gc(nsp)+bot
  283. GAMF=top/bot
  284. GM1=GAMF-1.0D0
  285. C-----------------------------------------
  286. C changing direction of the normal vector
  287. C-----------------------------------------
  288. CNX=-1.0D0*CNX
  289. CNY=-1.0D0*CNY
  290. IF(IDIM.EQ.2)THEN
  291. CTX=-1.0D0*CNY
  292. CTY=CNX
  293. ELSE
  294. CNZ=-1.0D0*CNZ
  295. CTX=-1.0D0*CTX
  296. CTY=-1.0D0*CTY
  297. CTZ=-1.0D0*CTZ
  298. CT2X=-1.0D0*CT2X
  299. CT2Y=-1.0D0*CT2Y
  300. CT2Z=-1.0D0*CT2Z
  301. ENDIF
  302. C---------------------------------------
  303. C******* On calcule UN, UT, UT2, ASON, S
  304. C---------------------------------------
  305. UNC=(UXC*CNX)+(UYC*CNY)+(UZC*CNZ)
  306. UNF=UNC
  307. UTC=(UXC*CTX)+(UYC*CTY)+(UZC*CTZ)
  308. C----------------------------------
  309. UXF=UNF*CNX+UTF*CTX+UT2F*CT2X
  310. UYF=UNF*CNY+UTF*CTY+UT2F*CT2Y
  311. UZF=UNF*CNZ+UTF*CTZ+UT2F*CT2Z
  312. C----------------------------------
  313. ECIN=0.5D0*((UXF*UXF)+(UYF*UYF)+(UZF*UZF))
  314. PSRF=(GM1/GAMF)*(HTF-ECIN)
  315. RF=PSRF/SF
  316. RF=RF**(1.0D0/GM1)
  317. PF=SF*(RF**GAMF)
  318. C-----------------------------------------------
  319. C******* Densite, vitesse, pression sur le bord
  320. C-----------------------------------------------
  321. MPRLI.VPOCHA(NLCB,1)=RF
  322. MPRLI.VPOCHA(NLCB,2)=UXF
  323. MPRLI.VPOCHA(NLCB,3)=UYF
  324. IF(IDIM.EQ.3) MPRLI.VPOCHA(NLCB,4)=UZF
  325. MPRLI.VPOCHA(NLCB,IDIM+2)=PF
  326. do 1040 i=1,(nsp-1)
  327. MPRLI.VPOCHA(NLCB,IDIM+2+I)=YF.YET(I)
  328. 1040 continue
  329. C---------------------------------------------------
  330. C******* Probleme de Riemann entre l'etat gauche
  331. C RF,UNC,UTF,UT2F,PF et l'etat droite
  332. C RC,UNC,UTC,UT2C,PC
  333. C On utilise AUSM+
  334. C Flux dans le repaire normale
  335. C---------------------------------------------------
  336. NESP=NSP-1
  337. IF(IDIM.EQ.2)THEN
  338. CALL FAUSMP(NESP,
  339. & GAMF,RF,PF,UNC,UTF,
  340. & GAMC,RC,PC,UNC,UTC,
  341. & YF.YET,YC.YET,
  342. & FLUX2D.FU,
  343. & CELLT)
  344. C-------------------------------------------------------
  345. C******* Residuum (son SPG a le meme ordre que MELEFC)
  346. C-------------------------------------------------------
  347. MPRES.VPOCHA(IFAC,1)=FLUX2D.FU(1)*SURF/VOLU
  348. MPRES.VPOCHA(IFAC,2)=((FLUX2D.FU(2)*CNX)+(FLUX2D.FU(3)*CTX))
  349. & *SURF/VOLU
  350. MPRES.VPOCHA(IFAC,3)=((FLUX2D.FU(2)*CNY)+(FLUX2D.FU(3)*CTY))
  351. & *SURF/VOLU
  352. MPRES.VPOCHA(IFAC,4)=FLUX2D.FU(4)*SURF/VOLU
  353. do 1050 i=1,(nsp-1)
  354. MPRES.VPOCHA(IFAC,4+I)=FLUX2D.FU(4+I)*SURF/VOLU
  355. 1050 continue
  356. ELSE
  357. CALL FAUSM3(NESP,
  358. & GAMF,RF,PF,UNC,UTF,UT2F,
  359. & GAMC,RC,PC,UNC,UTC,UT2C,
  360. & YF.YET,YC.YET,
  361. & FLUX3D.FU,
  362. & CELLT)
  363. C------------------------------------------------------
  364. C******* Residuum (son SPG a le meme ordre que MELEFC)
  365. C------------------------------------------------------
  366. MPRES.VPOCHA(IFAC,1)=FLUX3D.FU(1)*SURF/VOLU
  367. MPRES.VPOCHA(IFAC,2)=((FLUX3D.FU(2)*CNX)+(FLUX3D.FU(3)*CTX)+
  368. & (FLUX3D.FU(4)*CT2X))*SURF/VOLU
  369. MPRES.VPOCHA(IFAC,3)=((FLUX3D.FU(2)*CNY)+(FLUX3D.FU(3)*CTY)+
  370. & (FLUX3D.FU(4)*CT2Z))*SURF/VOLU
  371. MPRES.VPOCHA(IFAC,4)=((FLUX3D.FU(2)*CNZ)+(FLUX3D.FU(3)*CTZ)+
  372. & (FLUX3D.FU(4)*CT2Z))*SURF/VOLU
  373. MPRES.VPOCHA(IFAC,5)=FLUX3D.FU(5)*SURF/VOLU
  374. do 1060 i=1,(nsp-1)
  375. MPRES.VPOCHA(IFAC,5+I)=FLUX3D.FU(5+I)*SURF/VOLU
  376. 1060 continue
  377. ENDIF
  378. ENDIF
  379. 1 CONTINUE
  380. C
  381. SEGDES MELEFC
  382. C
  383. SEGDES MLEMC
  384. SEGDES MLEMCB
  385. SEGDES MLEMF
  386. C
  387. SEGDES MPNORM
  388. SEGDES MPVOL
  389. SEGDES MPSURF
  390. SEGDES MPRC
  391. SEGDES MPPC
  392. SEGDES MPVC
  393. SEGDES MPYN
  394. SEGDES MPLIM
  395. SEGDES MPRES
  396. SEGDES MPRLI
  397. SEGDES MLRECP
  398. SEGDES MLRECV
  399. SEGDES YC
  400. IF(TEMP .GT. 0.0D0) SEGDES YF
  401. SEGDES FLUX2D
  402. SEGDES FLUX3D
  403. C
  404. RETURN
  405. END
  406.  
  407.  
  408.  
  409.  
  410.  
  411.  
  412.  
  413.  
  414.  
  415.  
  416.  

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