1 : C FUNSHO SOURCE BECC 09/12/07 21:15:35 6579
2 : subroutine funsho(nordpo, Tmaxcv, acv_r, Rgas_r,
3 : & P_r, T_r, u_r, logr,
4 : & T, P, rho, u, d)
5 : C
6 : C INPUT
7 : C
8 : C nordpo = order of polynomial for cp and cv (see also
9 : C Tmaxcv)
10 : C
11 : C Tmaxcv = maximum temperature for cv polynomial expansion
12 : C cv(T) = cv(Tmaxcv) if T > Tmaxcv
13 : C
14 : C acv_r = to compute cv and ether for r (and rs);
15 : C vector such that
16 : C cv = \sum_{i=1,nordpo+1} acv(i) T^{i-1}
17 : C ether = \sum_{i=1,nordpo+1} acv(i) T^{i} / (i)
18 : C
19 : C Rgas_r = gas constant for r and (rs)
20 : C
21 : C p_r, T_r, u_r = pressure, temperature and velocities for r
22 : C
23 : C logr = .true. right travelling shock
24 : C = .false. left travelling shock
25 : C
26 : C T = temperature for rs
27 : C
28 : C OUTPUT
29 : C
30 : C p, u
31 : C rho, d = pressure, velociy, density and
32 : C shock speed for rs
33 : C NB. for a non-er shock, d = u + c
34 : C with c = speed sound in rs
35 : C
36 : c implicit none
37 : integer nordpo
38 : real*8 rindr
39 : real*8 Tmaxcv, acv_r(1:(nordpo+1)), Rgas_r
40 : & , P_r, T_r, u_r
41 : & , T, P, rho, u, d
42 : real*8 un_r, rho_r
43 : & , et, cv, et_r, cv_r, ht, ht_r, b, c, csi, b2
44 : logical logr
45 : C
46 : if (logr) then
47 : rindr = 1.0D0
48 : else
49 : rindr = -1.0D0
50 : endif
51 : un_r = rindr * u_r
52 : C
53 : rho_r = P_r / (Rgas_r * T_r)
53 : call prith1(nordpo, acv_r, Tmaxcv, T_r, et_r, cv_r)
54 : ht_r = et_r + (Rgas_r * T_r)
54 : call prith1(nordpo, acv_r, Tmaxcv, T, et, cv)
55 : ht = et + (Rgas_r * T)
56 : C
57 : b = ((ht - ht_r) / Rgas_r) + (0.5d0 * (T_r - T))
58 : c = T_r
59 : b2 = b * b
60 : if (((c * T) .lt. (1d-4 * b2)) .and. (b .lt. 0.0D0)) then
61 : csi = (0.5d0 * c / b2) - (0.125d0 * (T * ((c / b2)**2))) +
62 : & (0.0625d0 * (T*T * ((c / b2)**3)))
63 : csi = csi * abs(b)
64 : else
65 : csi = sqrt(b2 + (c * T)) + b
66 : csi = csi / T
67 : endif
68 : C
69 : rho = rho_R * csi
70 : P = Rgas_r * rho * T
71 : u = abs (2.0D0 * (ht - ht_r) *
72 : & (csi - 1.0D0) / (csi + 1.0D0))
73 : u = un_r + (sign(1.0D0,(csi - 1.0D0)) * sqrt(u))
74 : C
75 : if (T .lt. T_r) then
76 : d = (cv + Rgas_r) / cv
77 : d = d * Rgas_r * T
78 : d = u + sqrt(d)
79 : else
80 : if (abs (T - T_r) .lt. (1d-5 * T_r)) then
81 : d = 2.0D0 * (cv_r + Rgas_r) / cv_r
82 : d = d * Rgas_r * T_r * csi * csi / (csi + 1.0D0)
83 : else
84 : d = 2.0D0 * (ht - ht_r) * csi * csi / ((csi * csi) - 1.0D0)
85 : endif
86 : d = abs(d)
87 : d = un_r + sqrt(d)
88 : endif
89 : C
90 : u = rindr * u
91 : d = rindr * d
92 : C
93 : return
94 : end
95 :
© Cast3M 2003 - All rights reserved.
Disclaimer