/book/chap23/acoustics/rpt2acvq.f

http://github.com/clawpack/clawpack-4.x · FORTRAN Legacy · 95 lines · 50 code · 3 blank · 42 comment · 1 complexity · 82047e2add75e2c080427e6fea9f3490 MD5 · raw file

  1. c
  2. c
  3. c =====================================================
  4. subroutine rpt2(ixy,maxm,meqn,mwaves,mbc,mx,ql,qr,
  5. & aux1,aux2,aux3,imp,asdq,bmasdq,bpasdq)
  6. c =====================================================
  7. implicit double precision (a-h,o-z)
  8. c
  9. c # Riemann solver in the transverse direction for the 2d acoustics
  10. c # equations on a quadrilateral grid.
  11. c
  12. c # Split asdq (= A^* \Delta q, where * = + or -)
  13. c # into down-going flux difference bmasdqb (= B^- A^* \Delta q)
  14. c # and up-going flux difference bpasdq (= B^+ A^* \Delta q)
  15. c
  16. c
  17. dimension ql(1-mbc:maxm+mbc, meqn)
  18. dimension qr(1-mbc:maxm+mbc, meqn)
  19. dimension asdq(1-mbc:maxm+mbc, meqn)
  20. dimension bmasdq(1-mbc:maxm+mbc, meqn)
  21. dimension bpasdq(1-mbc:maxm+mbc, meqn)
  22. dimension aux1(1-mbc:maxm+mbc, 9)
  23. dimension aux2(1-mbc:maxm+mbc, 9)
  24. dimension aux3(1-mbc:maxm+mbc, 9)
  25. c
  26. c
  27. if (ixy.eq.1) then
  28. ilenrat = 3
  29. ixtran = 4
  30. iytran = 5
  31. else
  32. ilenrat = 6
  33. ixtran = 1
  34. iytran = 2
  35. endif
  36. c
  37. c
  38. do 20 i=2-mbc,mx+mbc
  39. c # impedance and sound speed in adjacent cells:
  40. zim = aux1(i,8)
  41. zi = aux2(i,8)
  42. zip = aux3(i,8)
  43. cim = aux1(i,9)
  44. ci = aux2(i,9)
  45. cip = aux3(i,9)
  46. c
  47. c # pressure component of asdq:
  48. c
  49. asdqp = asdq(i,1)
  50. c
  51. c
  52. c # up-going:
  53. c -----------
  54. c
  55. xtran = aux3(i,ixtran)
  56. ytran = aux3(i,iytran)
  57. c
  58. c # transverse velocity component of asdq:
  59. asdqt = xtran*asdq(i,2) + ytran*asdq(i,3)
  60. b2 = (asdqp + zi*asdqt) / (zi+zip)
  61. bpasdqp = b2*zip
  62. bpasdqt = b2
  63. c
  64. bpasdq(i,1) = cip * bpasdqp
  65. bpasdq(i,2) = cip * xtran*bpasdqt
  66. bpasdq(i,3) = cip * ytran*bpasdqt
  67. c
  68. c
  69. c # down-going:
  70. c -------------
  71. c
  72. xtran = aux2(i,ixtran)
  73. ytran = aux2(i,iytran)
  74. c
  75. c # transverse velocity component of asdq:
  76. asdqt = xtran*asdq(i,2) + ytran*asdq(i,3)
  77. b1 = (-asdqp + zi*asdqt) / (zim+zi)
  78. bmasdqp = -b1*zim
  79. bmasdqt = b1
  80. c
  81. bmasdq(i,1) = -cim * bmasdqp
  82. bmasdq(i,2) = -cim * xtran*bmasdqt
  83. bmasdq(i,3) = -cim * ytran*bmasdqt
  84. c
  85. c # multiply by ratio of lengths:
  86. do m=1,3
  87. bmasdq(i,m) = bmasdq(i,m)*aux2(i,9-ilenrat)
  88. bpasdq(i,m) = bpasdq(i,m)*aux3(i,9-ilenrat)
  89. enddo
  90. 20 continue
  91. return
  92. end