Télécharger funsho.eso

Retour à la liste

Numérotation des lignes :

  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)
  54. call prith1(nordpo, acv_r, Tmaxcv, T_r, et_r, cv_r)
  55. ht_r = et_r + (Rgas_r * T_r)
  56. call prith1(nordpo, acv_r, Tmaxcv, T, et, cv)
  57. ht = et + (Rgas_r * T)
  58. C
  59. b = ((ht - ht_r) / Rgas_r) + (0.5d0 * (T_r - T))
  60. c = T_r
  61. b2 = b * b
  62. if (((c * T) .lt. (1d-4 * b2)) .and. (b .lt. 0.0D0)) then
  63. csi = (0.5d0 * c / b2) - (0.125d0 * (T * ((c / b2)**2))) +
  64. & (0.0625d0 * (T*T * ((c / b2)**3)))
  65. csi = csi * abs(b)
  66. else
  67. csi = sqrt(b2 + (c * T)) + b
  68. csi = csi / T
  69. endif
  70. C
  71. rho = rho_R * csi
  72. P = Rgas_r * rho * T
  73. u = abs (2.0D0 * (ht - ht_r) *
  74. & (csi - 1.0D0) / (csi + 1.0D0))
  75. u = un_r + (sign(1.0D0,(csi - 1.0D0)) * sqrt(u))
  76. C
  77. if (T .lt. T_r) then
  78. d = (cv + Rgas_r) / cv
  79. d = d * Rgas_r * T
  80. d = u + sqrt(d)
  81. else
  82. if (abs (T - T_r) .lt. (1d-5 * T_r)) then
  83. d = 2.0D0 * (cv_r + Rgas_r) / cv_r
  84. d = d * Rgas_r * T_r * csi * csi / (csi + 1.0D0)
  85. else
  86. d = 2.0D0 * (ht - ht_r) * csi * csi / ((csi * csi) - 1.0D0)
  87. endif
  88. d = abs(d)
  89. d = un_r + sqrt(d)
  90. endif
  91. C
  92. u = rindr * u
  93. d = rindr * d
  94. C
  95. return
  96. end
  97.  
  98.  

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