Télécharger racre2.eso

Retour à la liste

Numérotation des lignes :

  1. C RACRE2 SOURCE BECC 09/12/07 21:16:00 6579
  2. subroutine racre2(nordpo, Tmaxcv,
  3. & Rgas_l, acv_l, Rgas_b, acv_b, Rgas_r, acv_r,
  4. & q, k0,
  5. & r_l, p_l, u_l, T_l, gam_l,
  6. & r_r, p_r, u_r, T_r, gam_r,
  7. & d_l, d_r,
  8. & r_d, p_d, u_d, T_d, d_d,
  9. & r_c, p_c, u_c, T_c, d_c,
  10. & r_b, p_b, u_b, T_b, d_b,
  11. & r_a, p_a, u_a, T_a, d_a,
  12. & logan, lognc)
  13. C
  14. c************************************************************************
  15. c
  16. c project :
  17. c
  18. c name : racre2
  19. c
  20. c description : computation of the shock-shock field-by-field
  21. c decomposition of the reactive Riemann problem for
  22. c thermally perfect gases.
  23. C
  24. C We have to distinguish between 6 different states
  25. C
  26. C l = left state;
  27. C burnt state (csi = csimax)
  28. C |
  29. C | Left GNL wave
  30. C |
  31. C d = left star state;
  32. C connected to l by a GNL wave
  33. C (csi = csimax)
  34. C |
  35. C | Contact discontinuity
  36. C |
  37. C c = right star star star state;
  38. C connected to rss state by a rarefaction
  39. C wave
  40. C (csi = csimax)
  41. C |
  42. C | Rarefaction wave (weak defl.)
  43. C |
  44. C b = right star star state;
  45. C connected to rs by a reactive wave
  46. C (csi = csimax)
  47. C |
  48. C | Reactive shock (weak defl.)
  49. C |
  50. C a = right star state;
  51. C connected to r by a GNL wave
  52. C (csi = csimin)
  53. C |
  54. C | Right GNL wave
  55. C |
  56. C r = right state;
  57. C unburnt state (csi = csimin)
  58. c
  59. c language : fortran 77
  60. c
  61. c author : a. beccantini den/dm2s/sfme/ltmf
  62. c
  63. c************************************************************************
  64. c
  65. c called by : fbecre
  66. c
  67. c called program : ...
  68. c
  69. c************************************************************************
  70. c
  71. c
  72. c**** input
  73. c
  74. c nordpo = polynomials degree of cv(t)
  75. c tmaxcv = threshold temperature for cv(t)
  76. c rgas_l = gas constant (J/kg/K)in l (and in ls)
  77. c acv_l = cv_l = \sum_{i=0,nordpo} (acv_l{i} T^i)
  78. c rgas_rss = gas constant in rss
  79. c acv_z = cv_rss = \sum_{i=0,nordpo} (acv_z{i} T^i)
  80. c rgas_r = gas constant r (and in rs)
  81. c acv_r = cv_r = \sum_{i=0,nordpo} (acv_r{i} T^i)
  82. c
  83. c q = released heat per mass unit from rs to rss
  84. c k0 = fundamental flame speed
  85. c
  86. C r_l = density in l
  87. C p_l = pressure in
  88. c u_l = velocity in l
  89. C T_l = temperature in l (redundant with r_l and
  90. c p_l)
  91. c gam_l = specific heat ratio in l (redundant
  92. c with r_l and p_l)
  93. c
  94. C r_r = density in r
  95. C p_r = pressure in r
  96. c u_r = velocity in r
  97. C T_r = temperature in r (redundant with r_r and
  98. c p_r)
  99. c gam_r = specific heat ratio in r (redundant
  100. c with r_r and p_r)
  101. c
  102. c**** output
  103. C
  104. C d_r, d_l = speed of states l and r
  105. c
  106. C r_d = density in d
  107. C p_d = pressure in d
  108. c u_d = velocity in d
  109. c T_d = temperature in d
  110. c d_d = shock speed of the left shock l<-d
  111. c
  112. C r_c = density in c
  113. C p_c = pressure in c
  114. c u_c = velocity in c
  115. c T_c = temperature in c
  116. c d_c = shock speed of the (non-entropy respecting)
  117. C right shock c->b
  118. c
  119. C r_b = density in b
  120. C p_b = pressure in b
  121. c u_b = velocity in b
  122. c T_b = temperature in b
  123. c d_b = shock speed of the reactive shock
  124. C right shock b->a
  125. c
  126. C r_a = density in a
  127. C p_a = pressure in a
  128. c u_a = velocity in a
  129. c T_a = temperature in a
  130. c d_a = shock speed of the
  131. C right shock a->l
  132. c
  133. c logan = an anomaly is detected
  134. c lognc = the approach has some problems
  135. C
  136. c
  137. c************************************************************************
  138. c
  139. c 07/12/2009 Created.
  140. c
  141. c implicit none
  142. integer nordpo
  143. integer nmax
  144. parameter (nmax = 20)
  145. integer iter, indreg
  146. real*8 eps
  147. parameter(eps = 1.0D-6)
  148. real*8 Tmaxcv, Rgas_l, Rgas_b, Rgas_r
  149. & , acv_l(0:nordpo), acv_b(0:nordpo), acv_r(0:nordpo)
  150. & , q, k0
  151. & , r_l, p_l, u_l, T_l, gam_l
  152. & , r_r, p_r, u_r, T_r, gam_r
  153. & , d_l, d_r
  154. & , r_d, p_d, u_d, T_d, d_d
  155. & , r_c, p_c, u_c, T_c, d_c
  156. & , r_b, p_b, u_b, T_b, d_b
  157. & , r_a, p_a, u_a, T_a, d_a
  158. real*8 k0cjdt, k0sdt
  159. & , d_acj, d_bcj
  160. & , r_acj, T_acj, p_acj, u_acj
  161. & , r_bcj, T_bcj, p_bcj, u_bcj
  162. & , P_asdt, r_asdt, u_asdt, d_sdt
  163. & , T_csdt, P_csdt, r_csdt, u_csdt
  164. & , T_cgp, P_cgp, r_cgp, u_cgp, d_cgp
  165. & , dt1, dt2
  166. C
  167. real*8 e_b, cv_b
  168. real*8 T_ag, T_bg, T_ag1, T_bg1, T_cg, T_dg
  169. real*8 Tref, pref, Dtmin, coefder
  170. real*8 P_ag, r_ag, u_ag, d_ag
  171. real*8 r_bg, p_bg, u_bg, d_bg
  172. real*8 r_cg, p_cg, u_cg, d_cg
  173. real*8 T_agp, P_agp, r_agp, u_agp, d_agp
  174. real*8 T_agpp, P_agpp, r_agpp, u_agpp, d_agpp
  175. real*8 r_bgp, p_bgp, u_bgp, T_bgp, d_bgp
  176. real*8 r_bgpp, p_bgpp, u_bgpp, T_bgpp, d_bgpp
  177. real*8 dtudtb
  178. real*8 dpduc, dpdud
  179. real*8 P_dg, r_dg, u_dg, d_dg
  180. real*8 T_dgp, P_dgp, r_dgp, u_dgp, d_dgp
  181. real*8 a, b
  182. logical lognc, logan, logdef, logde1, logdeb
  183. real*8 erro
  184. parameter(logdeb = .false.)
  185. CC
  186. CC**** We write the inputs
  187. CC
  188. C write(*,*) 'racre2.f'
  189. C write(*,'(I2,E12.4)') nordpo, tmaxcv
  190. C write(*,'(6E12.4)') rgas_l, (acv_l(icon), icon=0,nordpo)
  191. C write(*,'(6E12.4)') rgas_b, (acv_b(icon), icon=0,nordpo)
  192. C write(*,'(6E12.4)') rgas_r, (acv_r(icon), icon=0,nordpo)
  193. C write(*,'(2E12.4)') q, k0
  194. C write(*,'(6E12.4)') r_l, p_l, u_l, T_l, gam_l
  195. C write(*,'(6E12.4)') r_r, p_r, u_r, T_r, gam_r
  196. C stop
  197. C
  198. C**** Reference scales for a compressible solver.
  199. C For a low Mach number solver, the cut off speed is better than
  200. C the sound speed
  201. C
  202. Tref = 0.5d0 * ((T_l + T_r) + abs (T_l - T_r))
  203. Dtmin = 1.0D-3 * Tref
  204. coefder = 1.001D0
  205. pref = 0.5d0 * ((p_l + p_r) + abs (p_l - p_r))
  206. C
  207. C Initialisation of temperatures in a, b, c, d
  208. C
  209. T_a = T_r
  210. T_d = T_l
  211. C
  212. C T_b and T_c computed using AIbCC
  213. C i.e.
  214. C h_b = h_r + q
  215. c =>
  216. call prith1(nordpo, acv_b, Tmaxcv, T_l, e_b, cv_b)
  217. T_b = gam_r / (gam_r - 1.0D0) * Rgas_r * T_r
  218. T_b = T_b + q
  219. T_b = T_b / (cv_b + Rgas_b)
  220. T_c = T_b
  221. C
  222. C
  223. C TEST
  224. if (logdeb) then
  225. write(*,*) ' '
  226. write(*,*) 'racre2.f'
  227. write(*,*) 'Guess values'
  228. write(*,*) 'T_l, T_d, T_c, T_b, T_a, T_r '
  229. write(*,*) T_l, T_d, T_c, T_b, T_a, T_r
  230. endif
  231. C
  232. do iter = 1, nmax
  233. C
  234. C Guess temperatures at this iteration
  235. C
  236. T_ag = T_a
  237. T_bg = T_b
  238. T_cg = T_c
  239. T_dg = T_d
  240. C
  241. C states g1 to compute the properties (gamma)
  242. C
  243. T_ag1 = T_a
  244. T_bg1 = T_b
  245. C
  246. if (logdeb) then
  247. write(*,*) ' '
  248. write(*,*) 'iter = ', iter
  249. write(*,*) 'T_ag, T_bg, T_cg, T_dg ='
  250. write(*,*) T_ag, T_bg, T_cg, T_dg
  251. write(*,*) 'T_ag1, T_bg1 ='
  252. write(*,*) T_ag1, T_bg1
  253. endif
  254. C
  255. C******* Right state
  256. C
  257. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  258. CCCCCCCCC Test for funcj state CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  259. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  260. C acv computed for a mixture of h2-air at stoichiometric
  261. C condition (burnt vs unburn)
  262. C if (abs (Rgas_r - Rgas_b) .lt. (1.0D-6 * Rgas_b)) goto 9999
  263. C write(*,*) 'q =', q
  264. C p_r = 1.0D5
  265. C t_r = 290.0D0
  266. C u_r = 20.0D0
  267. C q = 3376364.59D0
  268. C k0 = 1000.0D0
  269. C T_ag1 = 1726.278895493752D0
  270. C T_ag = T_ag1
  271. C T_bg1 = 3519.399539515807D0
  272. C T_bg = T_bg1
  273. C T_cg = T_bg1 * 1.01D0
  274. CC
  275. C a) With these values and k0 = 1000 -> CJDT states
  276. C r_acj, T_acj, p_acj, u_acj, d_acj
  277. C 4.84002641182335 1665.81525651246
  278. C 3205701.36057718 1734.45181775933 2108.72324760979
  279. C r_bcj, T_bcj, p_bcj, u_bcj, d_bcj
  280. C 1.53143449865192 3355.53701100291 1740943.41687983
  281. C 925.856068333306 2108.72324760979
  282. C k0cjdt = 374.271429850455
  283. C
  284. C b) With these values and k0 = 100 -> CJDF states
  285. C r_acj, T_acj, p_acj, u_acj, d_acj
  286. C 1.89601907217341 420.523548481071 317016.244413484
  287. C 388.470810477263 699.102996003551
  288. C r_bcj, T_bcj, p_bcj, u_bcj, d_bcj
  289. C 0.193932466089758 2292.31572201026 150608.384363281
  290. C -489.198947662085 488.470810477263
  291. C k0cjdt = 374.271429850455
  292. C
  293. C c) With these values and k0 = 50 -> CJDF states
  294. C r_acj, T_acj, p_acj, u_acj, d_acj
  295. C 0.428757456680544 207.721167691767 35411.2237965742
  296. C -255.985416642147 85.8302953725479
  297. C r_bcj, T_bcj, p_bcj, u_bcj, d_bcj
  298. C 2.283207457276537E-002 2114.28097729292 16354.3117888704
  299. C -1144.92212016881 -205.985416642147
  300. C k0cjdt = 374.271429850455
  301. C
  302. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  303. CCCCCCCCC Test for funsdt state CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  304. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  305. C acv computed for a mixture of h2-air at stoichiometric
  306. C condition (burnt vs unburn)
  307. C if (abs (Rgas_r - Rgas_b) .lt. (1.0D-6 * Rgas_b)) goto 9999
  308. C write(*,*) 'q =', q
  309. C p_r = 1.0D5
  310. C t_r = 290.0D0
  311. C u_r = 20.0D0
  312. CC q = 3376364.59D0
  313. C k0 = 1000.0D0
  314. C T_ag1 = 1726.278895493752D0
  315. C T_ag = 1800.0D0
  316. C T_bg1 = 3519.399539515807D0
  317. C T_bg = T_bg1
  318. C T_cg = T_bg1 * 1.2D0
  319. C
  320. C With this value, for funcj we find the same values as for
  321. C case a.
  322. C Moreover, for funsdt, we find
  323. C T_r, P_r, u_r
  324. C 290.000000000000 100000.000000000 20.0000000000
  325. C T_ag, P_asdt, r_asdt, u_asdt, d_sdt, k0sdt
  326. C 1800.00000000000 3547944.85591171 4.95742167645132
  327. C 1831.11181823155 2215.13654202499 384.024723793441
  328. C T_csdt, P_csdt, r_csdt, u_csdt
  329. C 3625.49600977864 2494329.73747622 2.03077700146451
  330. C 1277.67640304670
  331. C
  332. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  333. CCCCCCCCC Test for funre1 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  334. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  335. CC
  336. CC acv computed for a mixture of h2-air at stoichiometric
  337. CC condition (burnt vs unburn)
  338. C if (abs (Rgas_r - Rgas_b) .lt. (1.0D-6 * Rgas_b)) goto 9999
  339. C write(*,*) 'q =', q
  340. C p_r = 1.0D5
  341. C t_r = 290.0D0
  342. C u_r = 20.0D0
  343. CCC q = 3376364.59D0
  344. C k0 = 100.0D0
  345. C T_ag1 = 1726.278895493752D0
  346. CC T_ag = 420.7D0
  347. C T_ag = 600.7D0
  348. C T_bg1 = 3519.399539515807D0
  349. C T_bg = T_bg1
  350. C T_cg = T_bg1 * 1.2D0
  351. CC
  352. CC With these values
  353. CC
  354. CC WDF
  355. CC T_r, P_r, u_r
  356. CC 290.000000000000 100000.000000000 20.0000000000000
  357. CC T_ag, P_ag, r_ag, u_ag, d_ag
  358. CC 600.700000000000 688354.959139395 2.88208190644972
  359. CC 708.663754615697 1005.09629332136
  360. CC T_bg, P_bg, r_bg, u_bg, d_bg
  361. CC 2669.18760214319 589828.308238861 0.652262358864745
  362. CC 366.804452560316 808.663754615697
  363. C
  364. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  365. CCCCCCCCC Test for fundt CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  366. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  367. CC
  368. CC acv computed for a mixture of h2-air at stoichiometric
  369. CC condition (burnt vs unburn)
  370. C if (abs (Rgas_r - Rgas_b) .lt. (1.0D-6 * Rgas_b)) goto 9999
  371. C write(*,*) 'q =', q
  372. C p_r = 1.0D5
  373. C t_r = 290.0D0
  374. C u_r = 20.0D0
  375. CC q = 3376364.59D0
  376. C k0 = 1000.0D0
  377. C T_ag1 = 1726.278895493752D0
  378. C T_ag = 600.7D0
  379. C T_bg1 = 3519.399539515807D0
  380. C T_bg = T_bg1
  381. C T_cg = 4000.0D0
  382. CC
  383. CC With these values
  384. CC
  385. CC SDT
  386. CC T_r, P_r, u_r
  387. CC 290.000000000000 100000.000000000 20.0000000000000
  388. CC T_ag, P_ag, r_ag, u_ag
  389. CC 2133.61520897134 4396903.39501966 5.18301338812827
  390. CC 2051.13020586146
  391. CC T_cg, P_cg, r_cg, u_cg, d_cg
  392. CC 4000.00000000000 3553808.96027808 2.62246614378061
  393. CC 1652.60261159855 2459.29495487122
  394. CC T_agp, P_agp, r_agp, u_agp
  395. CC 2137.48514727336 4406760.43821584 5.18522780966239
  396. CC 2053.54581926517
  397. CC T_cgp, P_cgp, r_cgp, u_cgp, d_cgp
  398. CC 4004.00000000000 3565189.78187338 2.62823616977283
  399. CC 1656.17695829123 2461.98642456893
  400. CC
  401. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  402. CCCCCCCCC END OF TEST CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  403. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  404. C
  405. C Computation of CJDF/CJDT states
  406. C
  407. call funcj(nordpo, Tmaxcv, acv_r, Rgas_r, acv_b, Rgas_b,
  408. & k0, q,
  409. & T_ag1, T_bg1,
  410. & p_r, T_r, u_r,
  411. & logdef,
  412. & r_acj, T_acj, p_acj, u_acj, d_acj,
  413. & r_bcj, T_bcj, p_bcj, u_bcj, d_bcj,
  414. & k0cjdt,
  415. & logan)
  416. C
  417. C Now I know the CJ state.
  418. C If logdef, CJ is the CJDF state.
  419. C Otherwise CJ is the CJDT state.
  420. C
  421. if (logan) then
  422. write(*,*) 'subroutine racre2.f'
  423. write(*,*) 'an error has been detected'
  424. goto 9999
  425. endif
  426. if (logdeb) then
  427. write(*,*) 'CJ state'
  428. write(*,*) 'Deflagration regime = ', logdef
  429. write(*,*) 'r_acj, T_acj, p_acj, u_acj, d_acj'
  430. write(*,*) r_acj, T_acj, p_acj, u_acj, d_acj
  431. write(*,*) 'r_bcj, T_bcj, p_bcj, u_bcj, d_bcj'
  432. write(*,*) r_bcj, T_bcj, p_bcj, u_bcj, d_bcj
  433. endif
  434. CC
  435. CC Test
  436. CC
  437. C write(*,*) 'deflagration ', logdef
  438. C write(*,*) ' r_acj, T_acj, p_acj, u_acj, d_acj '
  439. C write(*,*) r_acj, T_acj, p_acj, u_acj, d_acj
  440. C write(*,*) ' r_bcj, T_bcj, p_bcj, u_bcj, d_bcj '
  441. C write(*,*) r_bcj, T_bcj, p_bcj, u_bcj, d_bcj
  442. C write(*,*) 'k0cjdt = ', k0cjdt
  443. CC stop
  444. C
  445. if (T_cg .gt. T_bcj) then
  446. C
  447. C********** SDT/WDF
  448. C
  449. T_bg = max(T_bg, T_bcj)
  450. T_ag = max(T_ag, T_acj)
  451. C
  452. if (logdeb ) then
  453. write(*,*) 'WDF/SDT regime'
  454. write(*,*) 'We change T_ag, T_bg'
  455. write(*,*) T_ag, T_bg
  456. endif
  457. C
  458. logde1 = .true.
  459. C
  460. if( .not. logdef )then
  461. C We are here if k0 > k0cjdt
  462. C In this case T_ag >= T_a,cjdt
  463. C In this case T_bg >= T_b,cjdt
  464. C NB. We can have k0 > k0cjdt and being in WDF regime
  465. C
  466. C************* Computation of the SDT state (necessary to select the
  467. C regime) as function of T_ag
  468. C
  469. call funsdt(nordpo, Tmaxcv, acv_r, Rgas_r,
  470. & acv_b, Rgas_b, q, T_ag1, T_bg1,
  471. & P_r, T_r, u_r,
  472. & T_ag, P_asdt, r_asdt, u_asdt, d_sdt, k0sdt,
  473. & T_csdt, P_csdt, r_csdt, u_csdt,
  474. & logan)
  475. if (logan) then
  476. write(*,*) 'subroutine racre2.f'
  477. write(*,*) 'an error has been detected'
  478. goto 9999
  479. endif
  480. CC
  481. CC Test. We write the values thus determined
  482. CC
  483. C write(*,*) 'T_r, P_r, u_r'
  484. C write(*,*) T_r, P_r, u_r
  485. C write(*,*) 'T_ag, P_asdt, r_asdt, u_asdt, d_sdt, k0sdt'
  486. C write(*,*) T_ag, P_asdt, r_asdt, u_asdt, d_sdt, k0sdt
  487. C write(*,*) 'T_csdt, P_csdt, r_csdt, u_csdt'
  488. C write(*,*) T_csdt, P_csdt, r_csdt, u_csdt
  489. CC stop
  490. if (k0 .lt. k0sdt) then
  491. logde1 = .true.
  492. else
  493. logde1 = .false.
  494. endif
  495. endif
  496. C
  497. if (logde1) then
  498. indreg = 1
  499. C
  500. C************* WDF
  501. C
  502. C Computation of the state ag
  503. C
  504. call funsho( nordpo, Tmaxcv, acv_r, Rgas_r,
  505. & P_r, T_r, u_r, .true.,
  506. & T_ag, P_ag, r_ag, u_ag, d_ag)
  507. C
  508. C
  509. C Computation of the state bg (we change the value of T_bg)
  510. C
  511. call funre1(nordpo, Tmaxcv, Rgas_r, acv_r,
  512. & Rgas_b, acv_b, q, k0, T_ag1, T_bg1,
  513. & r_ag, p_ag, u_ag, T_ag,
  514. & r_bg, p_bg, u_bg, T_bg, d_bg, logan)
  515. if (logan) then
  516. write(*,*) 'subroutine racre2.f'
  517. write(*,*) 'an error has been detected'
  518. goto 9999
  519. endif
  520. C
  521. if (logdeb) then
  522. write(*,*) 'WDF'
  523. write(*,*) 'Parameter T_ag = ', T_ag
  524. write(*,*)
  525. & 'New values for T_bg = T_cg = ', T_bg
  526. endif
  527. C
  528. CC
  529. CC Test. We write the values thus determined
  530. CC
  531. C write(*,*) 'WDF'
  532. C write(*,*) 'T_r, P_r, u_r'
  533. C write(*,*) T_r, P_r, u_r
  534. C write(*,*) 'T_ag, P_ag, r_ag, u_ag, d_ag'
  535. C write(*,*) T_ag, P_ag, r_ag, u_ag, d_ag
  536. C write(*,*) 'T_bg, P_bg, r_bg, u_bg, d_bg'
  537. C write(*,*) T_bg, P_bg, r_bg, u_bg, d_bg
  538. C stop
  539. C
  540. if (d_bg .gt. d_ag) then
  541. if (logdeb) then
  542. write(*,*) 'WDF'
  543. write(*,*) 'd_bg, d_ag = ', d_bg, d_ag
  544. write(*,*) 'subroutine racre2.f'
  545. write(*,*) 'an error has been detected'
  546. endif
  547. C logan = .true.
  548. C goto 9999
  549. endif
  550. C
  551. C cg = bg
  552. C
  553. r_cg = r_bg
  554. p_cg = p_bg
  555. u_cg = u_bg
  556. T_cg = T_bg
  557. d_cg = d_bg
  558. C
  559. T_agp = max (T_ag + Dtmin, (T_ag * coefder))
  560. T_agpp = 2.0D0 * T_agp - T_ag
  561. C
  562. C Computation of states agp and agpp
  563. C
  564. call funsho( nordpo, Tmaxcv, acv_r, Rgas_r,
  565. & P_r, T_r, u_r, .true.,
  566. & T_agp, P_agp, r_agp, u_agp, d_agp)
  567. call funsho( nordpo, Tmaxcv, acv_r, Rgas_r,
  568. & P_r, T_r, u_r, .true.,
  569. & T_agpp, P_agpp, r_agpp, u_agpp, d_agpp)
  570. C
  571. call funre1(nordpo, Tmaxcv, Rgas_r, acv_r,
  572. & Rgas_b, acv_b, q, k0, T_ag1, T_bg1,
  573. & r_agp, p_agp, u_agp, T_agp,
  574. & r_bgp, p_bgp, u_bgp, T_bgp, d_bgp, logan)
  575. if (logan) then
  576. write(*,*) 'subroutine racre2.f'
  577. write(*,*) 'an error has been detected'
  578. goto 9999
  579. endif
  580. C
  581. call funre1(nordpo, Tmaxcv, Rgas_r, acv_r,
  582. & Rgas_b, acv_b, q, k0, T_ag1, T_bg1,
  583. & r_agpp, p_agpp, u_agpp, T_agpp,
  584. & r_bgpp, p_bgpp, u_bgpp, T_bgpp, d_bgpp, logan)
  585. if (logan) then
  586. write(*,*) 'subroutine racre2.f'
  587. write(*,*) 'an error has been detected'
  588. goto 9999
  589. endif
  590. C
  591. r_cgp = r_bgp
  592. p_cgp = p_bgp
  593. u_cgp = u_bgp
  594. T_cgp = T_bgp
  595. CC
  596. C write(*,*) 'WDF'
  597. C write(*,*) 'T_r, P_r, u_r'
  598. C write(*,*) T_r, P_r, u_r
  599. C write(*,*) 'T_agp, P_agp, r_agp, u_agp, d_agp'
  600. C write(*,*) T_agp, P_agp, r_agp, u_agp, d_agp
  601. C write(*,*) 'T_bgp, P_bgp, r_bgp, u_bgp, d_bgp'
  602. C write(*,*) T_bgp, P_bgp, r_bgp, u_bgp, d_bgp
  603. C write(*,*) 'T_agpp, P_agpp, r_agpp, u_agpp, d_agpp'
  604. C write(*,*) T_agpp, P_agpp, r_agpp, u_agpp, d_agpp
  605. C write(*,*) 'T_bgpp, P_bgpp, r_bgpp, u_bgpp, d_bgpp'
  606. C write(*,*) T_bgpp, P_bgpp, r_bgpp, u_bgpp, d_bgpp
  607. C
  608. C Second order formula (derivative evaluated in gp)
  609. C
  610. dTudTb = (T_agpp - T_ag) / (T_bgpp - T_bg)
  611. C
  612. C write(*,*) 'dTudTb = ', dTudTb
  613. C stop
  614. else
  615. indreg = 4
  616. C
  617. C************* SDT
  618. C Free parameter T_cg
  619. C Here we change the value of T_ag and T_bg
  620. C
  621. call fundt(nordpo, Tmaxcv, acv_r, Rgas_r,
  622. & acv_b, Rgas_b, q, T_ag1,
  623. & P_r, T_r, u_r,
  624. & T_cg, P_cg, r_cg, u_cg, d_cg, k0sdt,
  625. & T_ag, P_ag, r_ag, u_ag,
  626. & logan)
  627. C
  628. if (logan) then
  629. write(*,*) 'subroutine racre2.f'
  630. write(*,*) 'an error has been detected'
  631. goto 9999
  632. endif
  633. C
  634. if (logdeb) then
  635. write(*,*) 'SDT'
  636. write(*,*) 'Parameter T_cg = ', T_cg
  637. write(*,*) 'New values for T_ag = ', T_ag
  638. write(*,*) 'T_ag, P_ag, r_ag, u_ag'
  639. write(*,*) T_ag, P_ag, r_ag, u_ag
  640. write(*,*) 'T_cg, P_cg, r_cg, u_cg'
  641. write(*,*) T_cg, P_cg, r_cg, u_cg
  642. write(*,*) 'd_cg, k0sdt'
  643. write(*,*) d_cg, k0sdt
  644. endif
  645. C
  646. T_bg = T_cg
  647. P_bg = P_cg
  648. r_bg = r_cg
  649. u_bg = u_cg
  650. d_bg = d_cg
  651. d_ag = d_cg
  652. C
  653. T_cgp = max((T_cg + Dtmin), (T_cg * coefder))
  654. call fundt(nordpo, Tmaxcv, acv_r, Rgas_r,
  655. & acv_b, Rgas_b, q, T_ag1,
  656. & P_r, T_r, u_r,
  657. & T_cgp, P_cgp, r_cgp, u_cgp, d_cgp, k0sdt,
  658. & T_agp, P_agp, r_agp, u_agp,
  659. & logan)
  660. C
  661. if (logan) then
  662. write(*,*) 'subroutine racre2.f'
  663. write(*,*) 'an error has been detected'
  664. goto 9999
  665. endif
  666. C
  667. T_bgp = T_cgp
  668. P_bgp = P_cgp
  669. r_bgp = r_cgp
  670. u_bgp = u_cgp
  671. d_bgp = d_cgp
  672. d_agp = d_cgp
  673. CC
  674. CC Test. We write the values thus determined
  675. CC
  676. C write(*,*) 'SDT'
  677. C write(*,*) 'T_r, P_r, u_r'
  678. C write(*,*) T_r, P_r, u_r
  679. C write(*,*) 'T_ag, P_ag, r_ag, u_ag'
  680. C write(*,*) T_ag, P_ag, r_ag, u_ag
  681. C write(*,*) 'T_cg, P_cg, r_cg, u_cg, d_cg'
  682. C write(*,*) T_cg, P_cg, r_cg, u_cg, d_cg
  683. C write(*,*) 'T_agp, P_agp, r_agp, u_agp'
  684. C write(*,*) T_agp, P_agp, r_agp, u_agp
  685. C write(*,*) 'T_cgp, P_cgp, r_cgp, u_cgp, d_cgp'
  686. C write(*,*) T_cgp, P_cgp, r_cgp, u_cgp, d_cgp
  687. C stop
  688. C
  689. endif
  690. else
  691. C
  692. C********** CJDT/CJDF
  693. C
  694. r_ag = r_acj
  695. T_ag = T_acj
  696. P_ag = P_acj
  697. u_ag = u_acj
  698. d_ag = d_acj
  699. C
  700. r_bg = r_bcj
  701. T_bg = T_bcj
  702. P_bg = P_bcj
  703. u_bg = u_bcj
  704. d_bg = d_bcj
  705. C
  706. if (logdef) then
  707. C
  708. C CJDF
  709. C
  710. indreg = 2
  711. C
  712. if (logdeb) then
  713. write(*,*) 'CJDF'
  714. write(*,*) 'Parameter T_cg = ', T_cg
  715. write(*,*) 'New values for T_ag = ', T_ag
  716. endif
  717. C T_cg <= T_bg
  718. call funsho( nordpo, Tmaxcv, acv_b, Rgas_b,
  719. & P_bg, T_bg, u_bg, .true.,
  720. & T_cg, P_cg, r_cg, u_cg, d_cg)
  721. C
  722. C We take T_cgp > T_cg since a lower value gives problem
  723. C as T_cg goes to zero.
  724. C Note that as T_c changes, T_a and T_b do not
  725. C
  726. T_cgp = max((T_cg + Dtmin), (T_cg * coefder))
  727. call funsho( nordpo, Tmaxcv, acv_b, Rgas_b,
  728. & P_bg, T_bg, u_bg, .true.,
  729. & T_cgp, P_cgp, r_cgp, u_cgp, d_cgp)
  730. C
  731. C Even in this case, I need dTudTb and states
  732. C agp, agpp, bgp, bgpp
  733. C
  734. T_agp = max((T_ag + Dtmin), (T_ag*coefder))
  735. T_agpp = (2.0D0 * T_agp) - T_ag
  736. call funsho( nordpo, Tmaxcv, acv_b, Rgas_b,
  737. & P_r, T_r, u_r, .true.,
  738. & T_agp, P_agp, r_agp, u_agp, d_agp)
  739. call funsho( nordpo, Tmaxcv, acv_b, Rgas_b,
  740. & P_r, T_r, u_r, .true.,
  741. & T_agpp, P_agpp, r_agpp, u_agpp, d_agpp)
  742. C
  743. call funre1(nordpo, Tmaxcv, Rgas_r, acv_r,
  744. & Rgas_b, acv_b, q, k0, T_ag1, T_bg1,
  745. & r_agp, p_agp, u_agp, T_agp,
  746. & r_bgp, p_bgp, u_bgp, T_bgp, d_bgp, logan)
  747. if (logan) then
  748. write(*,*) 'subroutine racre2.f'
  749. write(*,*) 'an error has been detected'
  750. goto 9999
  751. endif
  752. call funre1(nordpo, Tmaxcv, Rgas_r, acv_r,
  753. & Rgas_b, acv_b, q, k0, T_ag1, T_bg1,
  754. & r_agpp, p_agpp, u_agpp, T_agpp,
  755. & r_bgpp, p_bgpp, u_bgpp, T_bgpp, d_bgpp, logan)
  756. if (logan) then
  757. write(*,*) 'subroutine racre2.f'
  758. write(*,*) 'an error has been detected'
  759. goto 9999
  760. endif
  761. C Second order formula (derivative evaluated in gp)
  762. dTudTb = (T_agpp - T_ag) / (T_bgpp - T_bg)
  763. else
  764. C
  765. C CJDT
  766. C
  767. indreg = 3
  768. C
  769. if (logdeb) then
  770. write(*,*) 'CJDT'
  771. write(*,*) 'Parameter T_cg = ', T_cg
  772. write(*,*) 'New values for T_ag = ', T_ag
  773. endif
  774. C T_cg <= T_bg
  775. call funsho( nordpo, Tmaxcv, acv_b, Rgas_b,
  776. & P_bg, T_bg, u_bg, .true.,
  777. & T_cg, P_cg, r_cg, u_cg, d_cg)
  778. C
  779. C We take T_cgp > T_cg.
  780. C Note that as T_c changes, T_a and T_b do not
  781. C
  782. T_cgp = max((T_cg + Dtmin), (T_cg * coefder))
  783. call funsho( nordpo, Tmaxcv, acv_b, Rgas_b,
  784. & P_bg, T_bg, u_bg, .true.,
  785. & T_cgp, P_cgp, r_cgp, u_cgp, d_cgp)
  786. C
  787. endif
  788. endif
  789. C
  790. dpduc = (P_cgp - P_cg) / (u_cgp - u_cg)
  791. C write(*,*) dpduc
  792. C
  793. C
  794. C******* Left state
  795. C
  796. call funsho( nordpo, Tmaxcv, acv_l, Rgas_l,
  797. & P_l, T_l, u_l, .false.,
  798. & T_dg, P_dg, r_dg, u_dg, d_dg)
  799. T_dgp = max((T_dg + Dtmin), (T_dg * coefder))
  800. call funsho( nordpo, Tmaxcv, acv_l, Rgas_l,
  801. & P_l, T_l, u_l, .false.,
  802. & T_dgp, P_dgp, r_dgp, u_dgp, d_dgp)
  803. if (dpduc .lt. 0.0d0) then
  804. write(*,*) ' dpduc < 0 , ', dpduc
  805. logan = .true.
  806. write(*,*) 'subroutine racre2.f'
  807. write(*,*) 'an error has been detected'
  808. goto 9999
  809. endif
  810. C
  811. dpdud = (P_dgp - P_dg) / (u_dgp - u_dg)
  812. if (dpdud .gt. 0.0d0) then
  813. write(*,*) 'dpdud > 0 , ', dpdud
  814. logan = .true.
  815. write(*,*) 'subroutine racre2.f'
  816. write(*,*) 'an error has been detected'
  817. goto 9999
  818. endif
  819. C
  820. C******* Intersection
  821. C
  822. a = dpduc - dpdud
  823. b = ((dpduc * u_cg) - (dpdud * u_dg)) + (p_dg - p_cg)
  824. u_d = b / a
  825. p_d = p_dg + (dpdud * (u_d - u_dg))
  826. p_c = p_cg + (dpduc * (u_d - u_cg))
  827. if ( abs (p_d - p_c) .gt. (eps * pref)) then
  828. write(*,*) 'pc - pd too big, ', (p_c - p_d)
  829. logan = .true.
  830. write(*,*) 'subroutine racre2.f'
  831. write(*,*) 'an error has been detected'
  832. goto 9999
  833. endif
  834. C p_c, p_d should be positive
  835. p_c = max (p_c, (eps * p_cg))
  836. p_d = max (p_d, (eps * p_dg))
  837. C
  838. C******* Computation of temperature T_d
  839. C
  840. T_d = T_dg + ((p_d - p_dg) * (T_dgp - T_dg) / (p_dgp - p_dg))
  841. T_d = max (T_d, (eps * T_dg))
  842. C
  843. T_c = T_cg + ((p_c - p_cg) * (T_cgp - T_cg) / (p_cgp - p_cg))
  844. T_c = max (T_c, (eps * T_cg))
  845. C
  846. if ((indreg .eq. 3) .or. (indreg .eq. 4)) then
  847. if( indreg .eq. 3) then
  848. C CJDT
  849. T_a = T_ag
  850. T_b = T_bg
  851. else
  852. T_b = T_c
  853. T_a = T_ag
  854. endif
  855. else
  856. C
  857. C indreg = 1 (WDF) or (indreg = 2) CJDF
  858. C
  859. if (T_c .le. T_bcj) then
  860. T_a = T_acj
  861. T_b = T_bcj
  862. else
  863. T_b = T_c
  864. b = dTudTb
  865. b = b / (2.0D0 * (T_bgp - T_bcj))
  866. a = -2.0D0 * b * (T_bcj - T_bg)
  867. dt1 = ((a + (b * (T_b - T_bg))) *
  868. & (T_b - T_bg))
  869. dt2 = dTudTb * (T_b - T_bg)
  870. if (abs(dt1) .lt. abs(dt2)) then
  871. T_a = T_ag + dt1
  872. else
  873. T_a = T_ag + dt2
  874. endif
  875. T_a = max (T_a, (eps * T_a))
  876. endif
  877. endif
  878. C
  879. C
  880. C******* TEST
  881. C
  882. erro = abs (T_a - T_ag) + abs (T_b - T_bg) +
  883. & abs(T_c - T_cg) + abs (T_d - T_dg) +
  884. & abs (T_a - T_ag1) + abs (T_b - T_bg1)
  885. erro = erro / Tref
  886. if (logdeb) then
  887. write(*,*) 'Intersection'
  888. write(*,*) 'T_ag1, T_bg1 = ', T_ag1, T_bg1
  889. write(*,*) 'dpduc, dpdud = ', dpduc, dpdud
  890. write(*,*) 'p_c, p_d = ', p_c, p_d
  891. write(*,*) 'T_d, T_c = ', T_d, T_c
  892. write(*,*) 'T_a, T_b = ', T_a, T_b
  893. WRITE(*,*) 'erro =', erro
  894. endif
  895. if (erro .Le. 1.0D-14) then
  896. lognc = .false.
  897. goto 9998
  898. endif
  899. enddo
  900. if (logdeb ) then
  901. write(*,*) 'Warning. Convergence non-reached'
  902. write(*,*) 'erro = ', erro
  903. endif
  904. lognc = .true.
  905. 9998 continue
  906. if (indreg .eq. 1) then
  907. call funsho( nordpo, Tmaxcv, acv_r, Rgas_r,
  908. & P_r, T_r, u_r, .true.,
  909. & T_a, P_a, r_a, u_a, d_a)
  910. call funre1(nordpo, Tmaxcv, Rgas_r, acv_r,
  911. & Rgas_b, acv_b, q, k0, T_a, T_b,
  912. & r_a, p_a, u_a, T_a,
  913. & r_b, p_b, u_b, T_b, d_b, logan)
  914. d_b = min(d_a,d_b)
  915. r_c = r_b
  916. p_c = p_b
  917. u_c = u_b
  918. T_c = T_b
  919. d_c = d_b
  920. endif
  921. if (indreg .eq. 4) then
  922. call fundt(nordpo, Tmaxcv, acv_r, Rgas_r,
  923. & acv_b, Rgas_b, q, T_a,
  924. & P_r, T_r, u_r,
  925. & T_c, P_c, r_c, u_c, d_c, k0sdt,
  926. & T_a, P_a, r_a, u_a,
  927. & logan)
  928. d_a = d_c
  929. r_b = r_c
  930. p_b = p_c
  931. u_b = u_c
  932. T_b = T_c
  933. d_b = d_c
  934. if (logdeb) then
  935. write(*,*) 'SDT'
  936. write(*,*) 'T_a, P_a, r_a, u_a'
  937. write(*,*) T_a, P_a, r_a, u_a
  938. write(*,*) 'T_c, P_c, r_c, u_c'
  939. write(*,*) T_c, P_c, r_c, u_c
  940. write(*,*) 'd_c, k0sdt'
  941. write(*,*) d_c, k0sdt
  942. endif
  943. endif
  944. if ((indreg .eq. 2) .or. (indreg .eq. 3)) then
  945. call funcj(nordpo, Tmaxcv, acv_r, Rgas_r, acv_b, Rgas_b,
  946. & k0, q,
  947. & T_a, T_b,
  948. & p_r, T_r, u_r,
  949. & logdef,
  950. & r_a, T_a, p_a, u_a, d_a,
  951. & r_b, T_b, p_b, u_b, d_b,
  952. & k0cjdt,
  953. & logan)
  954. if (logan) goto 9999
  955. call funsho( nordpo, Tmaxcv, acv_b, Rgas_b,
  956. & P_b, T_b, u_b, .true.,
  957. & T_c, P_c, r_c, u_c, d_c)
  958. d_b = min(d_a, d_b)
  959. d_c = min(d_b, d_c)
  960. endif
  961. call funsho( nordpo, Tmaxcv, acv_l, Rgas_l,
  962. & P_l, T_l, u_l, .false.,
  963. & T_d, P_d, r_d, u_d, d_d)
  964. d_d = min(d_c, d_d)
  965. d_l = gam_l * P_l / r_l
  966. d_l = u_l - sqrt(d_l)
  967. d_l = min(d_l, d_d)
  968. d_r = gam_r * P_r / r_r
  969. d_r = u_r + sqrt(d_r)
  970. d_r = max(d_r, d_a)
  971. C
  972. if (logdeb) then
  973. write(*,*) 'Final result'
  974. write(*,*) 'r'
  975. write(*,*) r_r, p_r, u_r, T_r, d_r
  976. write(*,*) 'a'
  977. write(*,*) r_a, p_a, u_a, T_a, d_a
  978. write(*,*) 'b'
  979. write(*,*) r_b, p_b, u_b, T_b, d_b
  980. write(*,*) 'c'
  981. write(*,*) r_c, p_c, u_c, T_c, d_c
  982. write(*,*) 'd'
  983. write(*,*) r_d, p_d, u_d, T_d, d_d
  984. write(*,*) 'l'
  985. write(*,*) r_l, p_l, u_l, T_l, d_l
  986. C write(*,*) 'acv'
  987. C do iter = 0, nordpo, 1
  988. C write(*,*) acv_b(iter), acv_l(iter)
  989. C enddo
  990. endif
  991. 9999 continue
  992. return
  993. end
  994.  
  995.  

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