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