DISPATCH
riemann_mod.f90
1 module riemann_mod
2  USE amr_parameters
3  USE hydro_parameters
4  USE const
5  USE io_mod
6  USE mpi_mod
7  USE trace_mod
8  implicit none
9  private
10  type riemann_t
11  integer(8):: ndiff=0, nsolu=0
12  real:: max_dlogd=10.
13  contains
14  procedure:: hlld
15  procedure, nopass:: hlld_2d
16  procedure, nopass:: llf
17  procedure, nopass:: llf_2d
18  end type
19  type(riemann_t), public:: rieman
20 contains
21 
22 !===============================================================================
23 !===============================================================================
24 subroutine hlld (self, qleft, qright, fgdnv, ngrid, detailed_timer)
25  class(riemann_t):: self
26  integer:: ngrid
27  real(dp),dimension(:,:)::qleft,qright,fgdnv
28  LOGICAL, OPTIONAL:: detailed_timer
29  REAL(dp)::sl,sr
30  REAL(dp)::entho,lim
31  REAL(dp)::rl,pl,ul,ecinl,etotl,ptotl,bl,cl,wstarl,vtotbl,emagl,sal,sqrrstarl,vl,wl
32  REAL(dp)::rr,pr,ur,ecinr,etotr,ptotr,br,cr,wstarr,vtotbr,emagr,sar,sqrrstarr,vr,wr
33  REAL(dp)::cfastl,rcl,rstarl,el,vdotbstarl,vstarl,etotstarstarl,etotstarl,vdotbl
34  REAL(dp)::cfastr,rcr,rstarr,er,vdotbstarr,vstarr,etotstarstarr,etotstarr,vdotbr
35  REAL(dp)::ustar,ptotstar,estar,vstarstar,wstarstar
36  REAL(dp)::ro,uo,ptoto,etoto,vo,wo
37  REAL(dp)::a,sgnm
38  REAL(dp)::bstarl,bstarr,cstarl,cstarr,calfvenl,calfvenr
39  REAL(dp)::bo,co,bstarstar,cstarstar,vdotbstarstar,vdotbo,emago
40 #define INLINE
41 #ifdef INLINE
42  REAL(dp)::d,p,u,v,w,b,c,b2,c2,d2,cf
43 #endif
44  INTEGER ::ivar,i,icase, ndiff
45  logical:: flag
46  integer, save:: itimer=0
47  !.............................................................................
48  !call trace%begin('riemann_t%hlld', itimer=itimer, detailed_timer=detailed_timer)
49  !
50  ! constants
51  if (isothermal) then
52  entho = one
53  else
54  entho = one/max(gamma-one,1e-6)
55  end if
56  !
57  fgdnv = 0d0
58  ndiff = 0
59  do i=1,ngrid
60  !print '(i3,1p,8e10.2,2x,8e10.2)',i,qleft(i,:),qright(i,:)
61  !
62  ! Enforce continuity of normal component
63  a=half*(qleft(i,4)+qright(i,4))
64  sgnm=sign(one,a)
65  qleft(i,4)=a; qright(i,4)=a
66  !
67  ! Left variables
68  rl=qleft(i,1); pl=qleft(i,2); ul=qleft(i,3)
69  vl=qleft(i,5); bl=qleft(i,6); wl=qleft(i,7); cl=qleft(i,8)
70  ecinl = half*(ul*ul+vl*vl+wl*wl)*rl
71  emagl = half*(a*a+bl*bl+cl*cl)
72  etotl = pl*entho+ecinl+emagl
73  ptotl = pl + emagl
74  vdotbl= ul*a+vl*bl+wl*cl
75  !
76  ! Right variables
77  rr=qright(i,1); pr=qright(i,2); ur=qright(i,3)
78  vr=qright(i,5); br=qright(i,6); wr=qright(i,7); cr=qright(i,8)
79  ecinr = half*(ur*ur+vr*vr+wr*wr)*rr
80  emagr = half*(a*a+br*br+cr*cr)
81  etotr = pr*entho+ecinr+emagr
82  ptotr = pr + emagr
83  vdotbr= ur*a+vr*br+wr*cr
84  !
85  ! Find the largest eigenvalues in the normal direction to the interface
86 #ifdef INLINE
87  d=qleft(i,1); p=qleft(i,2); u=qleft(i,3); a=qleft(i,4)
88  v=qleft(i,5); b=qleft(i,6); w=qleft(i,7); c=qleft(i,8)
89  b2 = a*a+b*b+c*c
90  c2 = gamma*p/d
91  d2 = half*(b2/d+c2)
92  cf = sqrt( d2 + sqrt(max(d2*d2-c2*a*a/d,zero)) )
93  cfastl=cf
94  d=qright(i,1); p=qright(i,2); u=qright(i,3); a=qright(i,4)
95  v=qright(i,5); b=qright(i,6); w=qright(i,7); c=qright(i,8)
96  b2 = a*a+b*b+c*c
97  c2 = gamma*p/d
98  d2 = half*(b2/d+c2)
99  cf = sqrt( d2 + sqrt(max(d2*d2-c2*a*a/d,zero)) )
100  cfastr=cf
101 #else
102  CALL find_speed_fast(qleft(i,:),cfastl)
103  CALL find_speed_fast(qright(i,:),cfastr)
104 #endif
105  !
106  ! Compute HLL wave speed
107  sl=min(ul,ur)-max(cfastl,cfastr)
108  sr=max(ul,ur)+max(cfastl,cfastr)
109  !
110  lim = 1.0_dp
111  !
112  ! Compute lagrangian sound speed
113  rcl=rl*(ul-sl)
114  rcr=rr*(sr-ur)
115  !
116  ! Compute acoustic star state
117  ustar =(rcr*ur +rcl*ul + (ptotl-ptotr))/(rcr+rcl)
118  ptotstar=(rcr*ptotl+rcl*ptotr+rcl*rcr*(ul-ur))/(rcr+rcl)
119  !
120  ! Left star region variables
121  rstarl=rl*(sl-ul)/(sl-ustar)
122  estar =rl*(sl-ul)*(sl-ustar)-a**2
123  el =rl*(sl-ul)*(sl-ul )-a**2
124  !if (dbg) print 1,'hlld4', estar,rl,SL,ul,ustar,A
125  if(abs(estar) < 1d-4*a**2)then
126  vstarl=vl
127  bstarl=bl
128  wstarl=wl
129  cstarl=cl
130  else
131  vstarl=vl-a*bl*(ustar-ul)/estar
132  bstarl=bl*el/estar
133  wstarl=wl-a*cl*(ustar-ul)/estar
134  cstarl=cl*el/estar
135  endif
136  vdotbstarl=ustar*a+vstarl*bstarl+wstarl*cstarl
137  etotstarl=((sl-ul)*etotl-ptotl*ul+ptotstar*ustar+a*(vdotbl-vdotbstarl))/(sl-ustar)
138  sqrrstarl=sqrt(rstarl)
139  calfvenl=abs(a)/sqrrstarl
140  sal=ustar-calfvenl
141  !
142  ! Right star region variables
143  rstarr=rr*(sr-ur)/(sr-ustar)
144  estar =rr*(sr-ur)*(sr-ustar)-a**2
145  er =rr*(sr-ur)*(sr-ur )-a**2
146  if(abs(estar) < 1d-4*a**2)then
147  vstarr=vr
148  bstarr=br
149  wstarr=wr
150  cstarr=cr
151  else
152  vstarr=vr-a*br*(ustar-ur)/estar
153  bstarr=br*er/estar
154  wstarr=wr-a*cr*(ustar-ur)/estar
155  cstarr=cr*er/estar
156  endif
157  !
158  vdotbstarr=ustar*a+vstarr*bstarr+wstarr*cstarr
159  etotstarr=((sr-ur)*etotr-ptotr*ur+ptotstar*ustar+a*(vdotbr-vdotbstarr))/(sr-ustar)
160  sqrrstarr=sqrt(rstarr)
161  calfvenr=abs(a)/sqrrstarr
162  sar=ustar+calfvenr
163  !
164  ! Double star region variables
165  vstarstar=(sqrrstarl*vstarl+sqrrstarr*vstarr+sgnm*(bstarr-bstarl))/(sqrrstarl+sqrrstarr)
166  wstarstar=(sqrrstarl*wstarl+sqrrstarr*wstarr+sgnm*(cstarr-cstarl))/(sqrrstarl+sqrrstarr)
167  bstarstar=(sqrrstarl*bstarr+sqrrstarr*bstarl+sgnm*sqrrstarl*sqrrstarr*(vstarr-vstarl))/(sqrrstarl+sqrrstarr)
168  cstarstar=(sqrrstarl*cstarr+sqrrstarr*cstarl+sgnm*sqrrstarl*sqrrstarr*(wstarr-wstarl))/(sqrrstarl+sqrrstarr)
169  vdotbstarstar=ustar*a+vstarstar*bstarstar+wstarstar*cstarstar
170  etotstarstarl=etotstarl-sgnm*sqrrstarl*(vdotbstarl-vdotbstarstar)
171  etotstarstarr=etotstarr+sgnm*sqrrstarr*(vdotbstarr-vdotbstarstar)
172 
173  ! Sample the solution at x/t=0
174  if(sl>0d0)then
175  icase=1
176  ro=rl
177  uo=ul
178  vo=vl
179  wo=wl
180  bo=bl
181  co=cl
182  ptoto=ptotl
183  etoto=etotl
184  vdotbo=vdotbl
185  emago=emagl
186  else if(sal>0d0)then
187  icase=2
188  ro=rstarl
189  uo=ustar
190  vo=vstarl
191  wo=wstarl
192  bo=bstarl
193  co=cstarl
194  ptoto=ptotstar
195  etoto=etotstarl
196  vdotbo=vdotbstarl
197  emago=half*(a*a+bo*bo+co*co)
198  else if(ustar>0d0)then
199  icase=3
200  ro=rstarl
201  uo=ustar
202  vo=vstarstar
203  wo=wstarstar
204  bo=bstarstar
205  co=cstarstar
206  ptoto=ptotstar
207  etoto=etotstarstarl
208  vdotbo=vdotbstarstar
209  emago=half*(a*a+bo*bo+co*co)
210  else if(sar>0d0)then
211  icase=4
212  ro=rstarr
213  uo=ustar
214  vo=vstarstar
215  wo=wstarstar
216  bo=bstarstar
217  co=cstarstar
218  ptoto=ptotstar
219  etoto=etotstarstarr
220  vdotbo=vdotbstarstar
221  emago=half*(a*a+bo*bo+co*co)
222  else if (sr>0d0)then
223  icase=5
224  ro=rstarr
225  uo=ustar
226  vo=vstarr
227  wo=wstarr
228  bo=bstarr
229  co=cstarr
230  ptoto=ptotstar
231  etoto=etotstarr
232  vdotbo=vdotbstarr
233  emago=half*(a*a+bo*bo+co*co)
234  else
235  icase=6
236  ro=rr
237  uo=ur
238  vo=vr
239  wo=wr
240  bo=br
241  co=cr
242  ptoto=ptotr
243  etoto=etotr
244  vdotbo=vdotbr
245  emago=emagr
246  end if
247  ! Add a diffusive flux in the down grad(ln(rho)) direction, which is only nonzero when the
248  ! smallest of the left/right density is within a factor 10 from smallr, and which has at most
249  ! magnitude of 0.25 times the largest of the left/right fast mode speeds.
250  ustar = alog(rl/rr)*max(0.0, alog10((10.*smallr)/min(rl,rr)))
251  ustar = 0.25*max(cfastr,cfastl)*min(1.0, max(-1.0, ustar))
252  uo = uo + ustar
253  ndiff = ndiff + merge(1,0,abs(ustar)>0.0)
254  !
255  ! Compute the Godunov flux
256  fgdnv(i,1) = ro*uo
257  fgdnv(i,2) = (etoto+ptoto)*uo-a*vdotbo
258  fgdnv(i,3) = ro*uo*uo+ptoto+emago*(lim-1.0_dp)-a*a*lim
259  fgdnv(i,4) = zero
260  fgdnv(i,5) = ro*uo*vo-a*bo*lim
261  fgdnv(i,6) = bo*uo-a*vo
262  fgdnv(i,7) = ro*uo*wo-a*co*lim
263  fgdnv(i,8) = co*uo-a*wo
264  end do
265  !$omp atomic
266  self%ndiff = self%ndiff + ndiff
267  !$omp atomic
268  self%nsolu = self%nsolu + ngrid
269  !call trace%end (itimer, detailed_timer=detailed_timer)
270 end subroutine hlld
271 
272 !===============================================================================
273 !> HLLD 2D fluxes
274 !===============================================================================
275 SUBROUTINE hlld_2d (qLL,qRL,qLR,qRR,emf,j,k,l0,l1,nvar,xdim,detailed_timer)
276  INTEGER, INTENT(IN) :: j,k,l0,l1,nvar,xdim
277  REAL(dp), DIMENSION(:,:),INTENT(IN)::qll,qrl,qlr,qrr
278  REAL(dp), DIMENSION(:,:,:):: emf
279  LOGICAL, OPTIONAL:: detailed_timer
280  ! local variables
281  INTEGER:: l
282  REAL(dp) :: e,ell,erl,elr,err,s,sl,sr,sb,st,sal,sar,sat,sab
283  REAL(dp) :: cfastllx,cfastrlx,cfastlrx,cfastrrx,cfastlly,cfastrly,cfastlry,cfastrry
284  REAL(dp) :: calfvenr,calfvenl,calfvent,calfvenb
285  REAL(dp) :: vllx,vrlx,vlrx,vrrx,vlly,vrly,vlry,vrry
286  REAL(dp) :: ptotll,ptotlr,ptotrl,ptotrr,rcllx,rclrx,rcrlx,rcrrx,rclly,rclry,rcrly,rcrry
287  REAL(dp) :: ustar,vstar,rstarllx,rstarlrx,rstarrlx,rstarrrx,rstarlly,rstarlry,rstarrly,rstarrry
288  REAL(dp) :: rstarll,rstarlr,rstarrl,rstarrr,astarll,astarlr,astarrl,astarrr,bstarll,bstarlr,bstarrl,bstarrr
289  REAL(dp) :: estarllx,estarlrx,estarrlx,estarrrx,estarlly,estarlry,estarrly,estarrry,estarll,estarlr,estarrl,estarrr
290  REAL(dp) :: astart,astarb,bstarr,bstarl
291 
292  REAL(dp) :: rll,rlr,rrl,rrr,pll,plr,prl,prr,ull,ulr,url,urr,vll,vlr,vrl,vrr, &
293  all,alr,arl,arr,bll,blr,brl,brr,cll,clr,crl,crr
294  REAL(dp), DIMENSION(1:nvar) :: qtmp
295  REAL(dp) :: b2, c2, d2, d1
296  !integer, save:: itimer(10)=0
297  integer, save:: itimer=0
298  !.............................................................................
299  !call trace%begin ('riemann_t%hlld_2d', itimer=itimer, detailed_timer=detailed_timer)
300  do l=l0,l1
301  !call trace%begin ('hlld_2d find_fast_speed', itimer=itimer(1))
302  rll=qll(l,1); pll=qll(l,2); ull=qll(l,3); vll=qll(l,4); all=qll(l,6); bll=qll(l,7) ; cll=qll(l,8)
303  rlr=qlr(l,1); plr=qlr(l,2); ulr=qlr(l,3); vlr=qlr(l,4); alr=qlr(l,6); blr=qlr(l,7) ; clr=qlr(l,8)
304  rrl=qrl(l,1); prl=qrl(l,2); url=qrl(l,3); vrl=qrl(l,4); arl=qrl(l,6); brl=qrl(l,7) ; crl=qrl(l,8)
305  rrr=qrr(l,1); prr=qrr(l,2); urr=qrr(l,3); vrr=qrr(l,4); arr=qrr(l,6); brr=qrr(l,7) ; crr=qrr(l,8)
306 
307  ! Compute 4 fast magnetosonic velocity relative to x direction
308  qtmp(1)=qll(l,1); qtmp(2)=qll(l,2); qtmp(7)=qll(l,5); qtmp(8)=qll(l,8)
309  qtmp(3)=qll(l,3); qtmp(4)=qll(l,6); qtmp(5)=qll(l,4); qtmp(6)=qll(l,7)
310  !call find_speed_fast(qtmp,cfastLLx)
311 
312  b2 = qtmp(4)**2+qtmp(6)**2+qtmp(8)**2
313  d1 = 1./qtmp(1)
314  c2 = gamma*qtmp(2)*d1
315  d2 = half*(b2*d1+c2)
316  cfastllx=sqrt(d2 + sqrt(max(d2*d2-c2*qtmp(4)**2*d1,zero)))
317 
318  qtmp(1)=qlr(l,1); qtmp(2)=qlr(l,2); qtmp(7)=qlr(l,5); qtmp(8)=qlr(l,8)
319  qtmp(3)=qlr(l,3); qtmp(4)=qlr(l,6); qtmp(5)=qlr(l,4); qtmp(6)=qlr(l,7)
320  !call find_speed_fast(qtmp,cfastLRx)
321 
322  b2 = qtmp(4)**2+qtmp(6)**2+qtmp(8)**2
323  d1 = 1./qtmp(1)
324  c2 = gamma*qtmp(2)*d1
325  d2 = half*(b2*d1+c2)
326  cfastlrx=sqrt(d2 + sqrt(max(d2*d2-c2*qtmp(4)**2*d1,zero)))
327 
328  qtmp(1)=qrl(l,1); qtmp(2)=qrl(l,2); qtmp(7)=qrl(l,5); qtmp(8)=qrl(l,8)
329  qtmp(3)=qrl(l,3); qtmp(4)=qrl(l,6); qtmp(5)=qrl(l,4); qtmp(6)=qrl(l,7)
330  !call find_speed_fast(qtmp,cfastRLx)
331 
332  b2 = qtmp(4)**2+qtmp(6)**2+qtmp(8)**2
333  d1 = 1./qtmp(1)
334  c2 = gamma*qtmp(2)*d1
335  d2 = half*(b2*d1+c2)
336 ! print*,'qtmp(4),qtmp(6),qtmp(8)',qtmp(4),qtmp(6),qtmp(8)
337  cfastrlx=sqrt(d2 + sqrt(max(d2*d2-c2*qtmp(4)**2*d1,zero)))
338 
339  qtmp(1)=qrr(l,1); qtmp(2)=qrr(l,2); qtmp(7)=qrr(l,5); qtmp(8)=qrr(l,8)
340  qtmp(3)=qrr(l,3); qtmp(4)=qrr(l,6); qtmp(5)=qrr(l,4); qtmp(6)=qrr(l,7)
341  qtmp(1)=max(qtmp(1),smallr)
342  !call find_speed_fast(qtmp,cfastRRx)
343 
344  b2 = qtmp(4)**2+qtmp(6)**2+qtmp(8)**2
345  d1 = 1./qtmp(1)
346  c2 = gamma*qtmp(2)*d1
347  d2 = half*(b2*d1+c2)
348  cfastrrx=sqrt(d2 + sqrt(max(d2*d2-c2*qtmp(4)**2*d1,zero)))
349 
350  ! Compute 4 fast magnetosonic velocity relative to y direction
351  qtmp(1)=qll(l,1); qtmp(2)=qll(l,2); qtmp(7)=qll(l,5); qtmp(8)=qll(l,8)
352  qtmp(3)=qll(l,4); qtmp(4)=qll(l,7); qtmp(5)=qll(l,3); qtmp(6)=qll(l,6)
353  !call find_speed_fast(qtmp,cfastLLy)
354 
355  b2 = qtmp(4)**2+qtmp(6)**2+qtmp(8)**2
356  d1 = 1./qtmp(1)
357  c2 = gamma*qtmp(2)*d1
358  d2 = half*(b2*d1+c2)
359  cfastlly=sqrt(d2 + sqrt(max(d2*d2-c2*qtmp(4)**2*d1,zero)))
360 
361  qtmp(1)=qlr(l,1); qtmp(2)=qlr(l,2); qtmp(7)=qlr(l,5); qtmp(8)=qlr(l,8)
362  qtmp(3)=qlr(l,4); qtmp(4)=qlr(l,7); qtmp(5)=qlr(l,3); qtmp(6)=qlr(l,6)
363  !call find_speed_fast(qtmp,cfastLRy)
364 
365  b2 = qtmp(4)**2+qtmp(6)**2+qtmp(8)**2
366  d1 = 1./qtmp(1)
367  c2 = gamma*qtmp(2)*d1
368  d2 = half*(b2*d1+c2)
369  cfastlry=sqrt(d2 + sqrt(max(d2*d2-c2*qtmp(4)**2*d1,zero)))
370 
371  qtmp(1)=qrl(l,1); qtmp(2)=qrl(l,2); qtmp(7)=qrl(l,5); qtmp(8)=qrl(l,8)
372  qtmp(3)=qrl(l,4); qtmp(4)=qrl(l,7); qtmp(5)=qrl(l,3); qtmp(6)=qrl(l,6)
373  !call find_speed_fast(qtmp,cfastRLy)
374 
375  b2 = qtmp(4)**2+qtmp(6)**2+qtmp(8)**2
376  d1 = 1./qtmp(1)
377  c2 = gamma*qtmp(2)*d1
378  d2 = half*(b2*d1+c2)
379  cfastrly=sqrt(d2 + sqrt(max(d2*d2-c2*qtmp(4)**2*d1,zero)))
380 
381  qtmp(1)=qrr(l,1); qtmp(2)=qrr(l,2); qtmp(7)=qrr(l,5); qtmp(8)=qrr(l,8)
382  qtmp(3)=qrr(l,4); qtmp(4)=qrr(l,7); qtmp(5)=qrr(l,3); qtmp(6)=qrr(l,6)
383  !call find_speed_fast(qtmp,cfastRRy)
384 
385  b2 = qtmp(4)**2+qtmp(6)**2+qtmp(8)**2
386  d1 = 1./qtmp(1)
387  c2 = gamma*qtmp(2)*d1
388  d2 = half*(b2*d1+c2)
389  cfastrry=sqrt(d2 + sqrt(max(d2*d2-c2*qtmp(4)**2*d1,zero)))
390 
391  sl=min(ull,ulr,url,urr)-max(cfastllx,cfastlrx,cfastrlx,cfastrrx)
392  sr=max(ull,ulr,url,urr)+max(cfastllx,cfastlrx,cfastrlx,cfastrrx)
393  sb=min(vll,vlr,vrl,vrr)-max(cfastlly,cfastlry,cfastrly,cfastrry)
394  st=max(vll,vlr,vrl,vrr)+max(cfastlly,cfastlry,cfastrly,cfastrry)
395  !call trace%end (itimer(1))
396 
397  if(sb>0d0)then
398  if(sl>0d0)then
399  !call trace%begin ('hlld_2d(2)', itimer=itimer(2))
400  ell=ull*bll-vll*all
401 
402  e=ell
403  !call trace%end (itimer(2))
404  else if(sr<0d0)then
405  !call trace%begin ('hlld_2d(3)', itimer=itimer(3))
406  erl=url*brl-vrl*arl
407 
408  e=erl
409  !call trace%end (itimer(3))
410  else
411  !call trace%begin ('hlld_2d(4)', itimer=itimer(4))
412  ptotll=pll+half*(all*all+bll*bll+cll*cll)
413  ptotlr=plr+half*(alr*alr+blr*blr+clr*clr)
414  ptotrl=prl+half*(arl*arl+brl*brl+crl*crl)
415  ptotrr=prr+half*(arr*arr+brr*brr+crr*crr)
416 
417  rcllx=rll*(ull-sl); rcrlx=rrl*(sr-url)
418  rclrx=rlr*(ulr-sl); rcrrx=rrr*(sr-urr)
419  rclly=rll*(vll-sb); rclry=rlr*(st-vlr)
420  rcrly=rrl*(vrl-sb); rcrry=rrr*(st-vrr)
421 
422  ustar=(rcllx*ull+rclrx*ulr+rcrlx*url+rcrrx*urr+(ptotll-ptotrl+ptotlr-ptotrr))/(rcllx+rclrx+rcrlx+rcrrx)
423  vstar=(rclly*vll+rclry*vlr+rcrly*vrl+rcrry*vrr+(ptotll-ptotlr+ptotrl-ptotrr))/(rclly+rclry+rcrly+rcrry)
424 
425  rstarllx=rll*(sl-ull)/(sl-ustar)
426  rstarll =rll*(sl-ull)/(sl-ustar)*(sb-vll)/(sb-vstar)
427  astarll=all*(sb-vll)/(sb-vstar)
428 
429  rstarlrx=rlr*(sl-ulr)/(sl-ustar)
430  rstarlr =rlr*(sl-ulr)/(sl-ustar)*(st-vlr)/(st-vstar)
431  astarlr=alr*(st-vlr)/(st-vstar)
432 
433  rstarrlx=rrl*(sr-url)/(sr-ustar)
434  rstarrl =rrl*(sr-url)/(sr-ustar)*(sb-vrl)/(sb-vstar)
435  astarrl=arl*(sb-vrl)/(sb-vstar)
436 
437  rstarrrx=rrr*(sr-urr)/(sr-ustar)
438  rstarrr =rrr*(sr-urr)/(sr-ustar)*(st-vrr)/(st-vstar)
439  astarrr=arr*(st-vrr)/(st-vstar)
440 
441  calfvenl=max(abs(alr)/sqrt(rstarlrx),abs(astarlr)/sqrt(rstarlr), &
442  abs(all)/sqrt(rstarllx),abs(astarll)/sqrt(rstarll),smallc)
443  calfvenr=max(abs(arr)/sqrt(rstarrrx),abs(astarrr)/sqrt(rstarrr), &
444  abs(arl)/sqrt(rstarrlx),abs(astarrl)/sqrt(rstarrl),smallc)
445 
446  sal=min(ustar-calfvenl,zero)
447  sar=max(ustar+calfvenr,zero)
448 
449  bstarll=bll*(sl-ull)/(sl-ustar)
450  bstarrl=brl*(sr-url)/(sr-ustar)
451 
452  estarllx=ustar*bstarll-vll*all
453  estarrlx=ustar*bstarrl-vrl*arl
454 
455  e=(sar*estarllx-sal*estarrlx+sar*sal*(brl-bll))/(sar-sal)
456  !call trace%end (itimer(4))
457  endif
458  else if (st<0d0)then
459  if(sl>0d0)then
460  !call trace%begin ('hlld_2d(5)', itimer=itimer(5))
461  elr=ulr*blr-vlr*alr
462 
463  e=elr
464  !call trace%end (itimer(5))
465  else if(sr<0d0)then
466  !call trace%begin ('hlld_2d(6)', itimer=itimer(6))
467  err=urr*brr-vrr*arr
468 
469  e=err
470  !call trace%end (itimer(6))
471  else
472  !call trace%begin ('hlld_2d(7)', itimer=itimer(7))
473  ptotll=pll+half*(all*all+bll*bll+cll*cll)
474  ptotlr=plr+half*(alr*alr+blr*blr+clr*clr)
475  ptotrl=prl+half*(arl*arl+brl*brl+crl*crl)
476  ptotrr=prr+half*(arr*arr+brr*brr+crr*crr)
477 
478  rcllx=rll*(ull-sl); rcrlx=rrl*(sr-url)
479  rclrx=rlr*(ulr-sl); rcrrx=rrr*(sr-urr)
480  rclly=rll*(vll-sb); rclry=rlr*(st-vlr)
481  rcrly=rrl*(vrl-sb); rcrry=rrr*(st-vrr)
482 
483  ustar=(rcllx*ull+rclrx*ulr+rcrlx*url+rcrrx*urr+(ptotll-ptotrl+ptotlr-ptotrr))/(rcllx+rclrx+rcrlx+rcrrx)
484  vstar=(rclly*vll+rclry*vlr+rcrly*vrl+rcrry*vrr+(ptotll-ptotlr+ptotrl-ptotrr))/(rclly+rclry+rcrly+rcrry)
485 
486  rstarllx=rll*(sl-ull)/(sl-ustar)
487  rstarll =rll*(sl-ull)/(sl-ustar)*(sb-vll)/(sb-vstar)
488  astarll=all*(sb-vll)/(sb-vstar)
489 
490  rstarlrx=rlr*(sl-ulr)/(sl-ustar)
491  rstarlr =rlr*(sl-ulr)/(sl-ustar)*(st-vlr)/(st-vstar)
492  astarlr=alr*(st-vlr)/(st-vstar)
493 
494  rstarrlx=rrl*(sr-url)/(sr-ustar)
495  rstarrl =rrl*(sr-url)/(sr-ustar)*(sb-vrl)/(sb-vstar)
496  astarrl=arl*(sb-vrl)/(sb-vstar)
497 
498  rstarrrx=rrr*(sr-urr)/(sr-ustar)
499  rstarrr =rrr*(sr-urr)/(sr-ustar)*(st-vrr)/(st-vstar)
500  astarrr=arr*(st-vrr)/(st-vstar)
501 
502  calfvenl=max(abs(alr)/sqrt(rstarlrx),abs(astarlr)/sqrt(rstarlr), &
503  abs(all)/sqrt(rstarllx),abs(astarll)/sqrt(rstarll),smallc)
504  calfvenr=max(abs(arr)/sqrt(rstarrrx),abs(astarrr)/sqrt(rstarrr), &
505  abs(arl)/sqrt(rstarrlx),abs(astarrl)/sqrt(rstarrl),smallc)
506 
507  sal=min(ustar-calfvenl,zero)
508  sar=max(ustar+calfvenr,zero)
509 
510  bstarlr=blr*(sl-ulr)/(sl-ustar)
511  bstarrr=brr*(sr-urr)/(sr-ustar)
512 
513  estarlrx=ustar*bstarlr-vlr*alr
514  estarrrx=ustar*bstarrr-vrr*arr
515 
516  e=(sar*estarlrx-sal*estarrrx+sar*sal*(brr-blr))/(sar-sal)
517  !call trace%end (itimer(7))
518  endif
519  else if(sl>0d0)then
520  !call trace%begin ('hlld_2d(8)', itimer=itimer(8))
521  ptotll=pll+half*(all*all+bll*bll+cll*cll)
522  ptotlr=plr+half*(alr*alr+blr*blr+clr*clr)
523  ptotrl=prl+half*(arl*arl+brl*brl+crl*crl)
524  ptotrr=prr+half*(arr*arr+brr*brr+crr*crr)
525 
526  rcllx=rll*(ull-sl); rcrlx=rrl*(sr-url)
527  rclrx=rlr*(ulr-sl); rcrrx=rrr*(sr-urr)
528  rclly=rll*(vll-sb); rclry=rlr*(st-vlr)
529  rcrly=rrl*(vrl-sb); rcrry=rrr*(st-vrr)
530 
531  ustar=(rcllx*ull+rclrx*ulr+rcrlx*url+rcrrx*urr+(ptotll-ptotrl+ptotlr-ptotrr))/(rcllx+rclrx+rcrlx+rcrrx)
532  vstar=(rclly*vll+rclry*vlr+rcrly*vrl+rcrry*vrr+(ptotll-ptotlr+ptotrl-ptotrr))/(rclly+rclry+rcrly+rcrry)
533 
534  rstarlly=rll*(sb-vll)/(sb-vstar)
535  rstarll =rll*(sl-ull)/(sl-ustar)*(sb-vll)/(sb-vstar)
536  bstarll=bll*(sl-ull)/(sl-ustar)
537 
538  rstarlry=rlr*(st-vlr)/(st-vstar)
539  rstarlr =rlr*(sl-ulr)/(sl-ustar)*(st-vlr)/(st-vstar)
540  bstarlr=blr*(sl-ulr)/(sl-ustar)
541 
542  rstarrly=rrl*(sb-vrl)/(sb-vstar)
543  rstarrl =rrl*(sr-url)/(sr-ustar)*(sb-vrl)/(sb-vstar)
544  bstarrl=brl*(sr-url)/(sr-ustar)
545 
546  rstarrry=rrr*(st-vrr)/(st-vstar)
547  rstarrr =rrr*(sr-urr)/(sr-ustar)*(st-vrr)/(st-vstar)
548  bstarrr=brr*(sr-urr)/(sr-ustar)
549 
550  calfvenb=max(abs(bll)/sqrt(rstarlly),abs(bstarll)/sqrt(rstarll), &
551  abs(brl)/sqrt(rstarrly),abs(bstarrl)/sqrt(rstarrl),smallc)
552  calfvent=max(abs(blr)/sqrt(rstarlry),abs(bstarlr)/sqrt(rstarlr), &
553  abs(brr)/sqrt(rstarrry),abs(bstarrr)/sqrt(rstarrr),smallc)
554 
555  sab=min(vstar-calfvenb,zero)
556  sat=max(vstar+calfvent,zero)
557 
558  astarll=all*(sb-vll)/(sb-vstar)
559  astarlr=alr*(st-vlr)/(st-vstar)
560 
561  estarlly=ull*bll-vstar*astarll
562  estarlry=ulr*blr-vstar*astarlr
563 
564  e=(sat*estarlly-sab*estarlry-sat*sab*(alr-all))/(sat-sab)
565  !call trace%end (itimer(8))
566  else if(sr<0d0)then
567  !call trace%begin ('hlld_2d(9)', itimer=itimer(9))
568  ptotll=pll+half*(all*all+bll*bll+cll*cll)
569  ptotlr=plr+half*(alr*alr+blr*blr+clr*clr)
570  ptotrl=prl+half*(arl*arl+brl*brl+crl*crl)
571  ptotrr=prr+half*(arr*arr+brr*brr+crr*crr)
572 
573  rcllx=rll*(ull-sl); rcrlx=rrl*(sr-url)
574  rclrx=rlr*(ulr-sl); rcrrx=rrr*(sr-urr)
575  rclly=rll*(vll-sb); rclry=rlr*(st-vlr)
576  rcrly=rrl*(vrl-sb); rcrry=rrr*(st-vrr)
577 
578  ustar=(rcllx*ull+rclrx*ulr+rcrlx*url+rcrrx*urr+(ptotll-ptotrl+ptotlr-ptotrr))/(rcllx+rclrx+rcrlx+rcrrx)
579  vstar=(rclly*vll+rclry*vlr+rcrly*vrl+rcrry*vrr+(ptotll-ptotlr+ptotrl-ptotrr))/(rclly+rclry+rcrly+rcrry)
580 
581  rstarlly=rll*(sb-vll)/(sb-vstar)
582  rstarll =rll*(sl-ull)/(sl-ustar)*(sb-vll)/(sb-vstar)
583  bstarll=bll*(sl-ull)/(sl-ustar)
584 
585  rstarlry=rlr*(st-vlr)/(st-vstar)
586  rstarlr =rlr*(sl-ulr)/(sl-ustar)*(st-vlr)/(st-vstar)
587  bstarlr=blr*(sl-ulr)/(sl-ustar)
588 
589  rstarrly=rrl*(sb-vrl)/(sb-vstar)
590  rstarrl =rrl*(sr-url)/(sr-ustar)*(sb-vrl)/(sb-vstar)
591  bstarrl=brl*(sr-url)/(sr-ustar)
592 
593  rstarrry=rrr*(st-vrr)/(st-vstar)
594  rstarrr =rrr*(sr-urr)/(sr-ustar)*(st-vrr)/(st-vstar)
595  bstarrr=brr*(sr-urr)/(sr-ustar)
596 
597  calfvenb=max(abs(bll)/sqrt(rstarlly),abs(bstarll)/sqrt(rstarll), &
598  abs(brl)/sqrt(rstarrly),abs(bstarrl)/sqrt(rstarrl),smallc)
599  calfvent=max(abs(blr)/sqrt(rstarlry),abs(bstarlr)/sqrt(rstarlr), &
600  abs(brr)/sqrt(rstarrry),abs(bstarrr)/sqrt(rstarrr),smallc)
601 
602  sab=min(vstar-calfvenb,zero)
603  sat=max(vstar+calfvent,zero)
604 
605  astarrl=arl*(sb-vrl)/(sb-vstar)
606  astarrr=arr*(st-vrr)/(st-vstar)
607 
608  estarrly=url*brl-vstar*astarrl
609  estarrry=urr*brr-vstar*astarrr
610 
611  e=(sat*estarrly-sab*estarrry-sat*sab*(arr-arl))/(sat-sab)
612  !call trace%end (itimer(9))
613  else
614  !call trace%begin ('hlld_2d(10)', itimer=itimer(10))
615  ptotll=pll+half*(all*all+bll*bll+cll*cll)
616  ptotlr=plr+half*(alr*alr+blr*blr+clr*clr)
617  ptotrl=prl+half*(arl*arl+brl*brl+crl*crl)
618  ptotrr=prr+half*(arr*arr+brr*brr+crr*crr)
619 
620  rcllx=rll*(ull-sl); rcrlx=rrl*(sr-url)
621  rclrx=rlr*(ulr-sl); rcrrx=rrr*(sr-urr)
622  rclly=rll*(vll-sb); rclry=rlr*(st-vlr)
623  rcrly=rrl*(vrl-sb); rcrry=rrr*(st-vrr)
624 
625  ustar=(rcllx*ull+rclrx*ulr+rcrlx*url+rcrrx*urr+(ptotll-ptotrl+ptotlr-ptotrr))/(rcllx+rclrx+rcrlx+rcrrx)
626  vstar=(rclly*vll+rclry*vlr+rcrly*vrl+rcrry*vrr+(ptotll-ptotlr+ptotrl-ptotrr))/(rclly+rclry+rcrly+rcrry)
627 
628  rstarllx=rll*(sl-ull)/(sl-ustar); bstarll=bll*(sl-ull)/(sl-ustar)
629  rstarlly=rll*(sb-vll)/(sb-vstar); astarll=all*(sb-vll)/(sb-vstar)
630  rstarll =rll*(sl-ull)/(sl-ustar)*(sb-vll)/(sb-vstar)
631 
632  rstarlrx=rlr*(sl-ulr)/(sl-ustar); bstarlr=blr*(sl-ulr)/(sl-ustar)
633  rstarlry=rlr*(st-vlr)/(st-vstar); astarlr=alr*(st-vlr)/(st-vstar)
634  rstarlr =rlr*(sl-ulr)/(sl-ustar)*(st-vlr)/(st-vstar)
635 
636  rstarrlx=rrl*(sr-url)/(sr-ustar); bstarrl=brl*(sr-url)/(sr-ustar)
637  rstarrly=rrl*(sb-vrl)/(sb-vstar); astarrl=arl*(sb-vrl)/(sb-vstar)
638  rstarrl =rrl*(sr-url)/(sr-ustar)*(sb-vrl)/(sb-vstar)
639 
640  rstarrrx=rrr*(sr-urr)/(sr-ustar); bstarrr=brr*(sr-urr)/(sr-ustar)
641  rstarrry=rrr*(st-vrr)/(st-vstar); astarrr=arr*(st-vrr)/(st-vstar)
642  rstarrr =rrr*(sr-urr)/(sr-ustar)*(st-vrr)/(st-vstar)
643 
644  calfvenl=max(abs(alr)/sqrt(rstarlrx),abs(astarlr)/sqrt(rstarlr), &
645  abs(all)/sqrt(rstarllx),abs(astarll)/sqrt(rstarll),smallc)
646  calfvenr=max(abs(arr)/sqrt(rstarrrx),abs(astarrr)/sqrt(rstarrr), &
647  abs(arl)/sqrt(rstarrlx),abs(astarrl)/sqrt(rstarrl),smallc)
648  calfvenb=max(abs(bll)/sqrt(rstarlly),abs(bstarll)/sqrt(rstarll), &
649  abs(brl)/sqrt(rstarrly),abs(bstarrl)/sqrt(rstarrl),smallc)
650  calfvent=max(abs(blr)/sqrt(rstarlry),abs(bstarlr)/sqrt(rstarlr), &
651  abs(brr)/sqrt(rstarrry),abs(bstarrr)/sqrt(rstarrr),smallc)
652  sal=min(ustar-calfvenl,zero)
653  sar=max(ustar+calfvenr,zero)
654  sab=min(vstar-calfvenb,zero)
655  sat=max(vstar+calfvent,zero)
656  estarll =ustar*bstarll-vstar*astarll
657  estarrl =ustar*bstarrl-vstar*astarrl
658  estarlr =ustar*bstarlr-vstar*astarlr
659  estarrr =ustar*bstarrr-vstar*astarrr
660 
661  astart=(sar*astarrr-sal*astarlr)/(sar-sal); astarb=(sar*astarrl-sal*astarll)/(sar-sal)
662  bstarr=(sat*bstarrr-sab*bstarrl)/(sat-sab); bstarl=(sat*bstarlr-sab*bstarll)/(sat-sab)
663 
664  e=(sal*sab*estarrr-sal*sat*estarrl-sar*sab*estarlr+sar*sat*estarll)/(sar-sal)/(sat-sab) &
665  & -sat*sab/(sat-sab)*(astart-astarb)+sar*sal/(sar-sal)*(bstarr-bstarl)
666  !call trace%end (itimer(10))
667  endif
668  if (xdim == 1) then
669  emf(l,j,k) = e
670  else if (xdim == 2) then
671  emf(k,l,j) = e
672  else if (xdim == 3) then
673  emf(j,k,l) = e
674  endif
675  end do
676  !call trace%end (itimer, detailed_timer=detailed_timer)
677 END SUBROUTINE hlld_2d
678 
679 !===============================================================================
680 !> LLF 2D fluxes
681 !===============================================================================
682 SUBROUTINE llf_2d (qLL,qRL,qLR,qRR,emf,j,k,l0,l1,nvar,xdim,detailed_timer)
683  INTEGER, INTENT(IN) :: j,k,l0,l1,nvar,xdim
684  REAL(dp), DIMENSION(:,:),INTENT(IN)::qll,qrl,qlr,qrr
685  REAL(dp), DIMENSION(:,:,:):: emf
686  LOGICAL, OPTIONAL:: detailed_timer
687  ! local variables
688  INTEGER:: l
689  REAL(dp) :: e, ev, ell, erl, elr, err, qleft(8), qright(8), fmean_x(8), fmean_y(8)
690  integer, save:: itimer=0
691  !.............................................................................
692  !call trace%begin ('riemann_t%hlld_2d', itimer=itimer, detailed_timer=detailed_timer)
693  !
694  do l=l0,l1
695  ! vx*by - vy*bx at the four edge centers
696  ell = qll(l,3)*qll(l,7) - qll(l,4)*qll(l,6)
697  erl = qrl(l,3)*qrl(l,7) - qrl(l,4)*qrl(l,6)
698  elr = qlr(l,3)*qlr(l,7) - qlr(l,4)*qlr(l,6)
699  err = qrr(l,3)*qrr(l,7) - qrr(l,4)*qrr(l,6)
700  !
701  ! find the average value of E
702  ev = forth*(ell+erl+elr+err)
703  !
704  ! call the first solver in the x direction
705  ! density
706  qleft(1) = half*(qll(l,1)+qlr(l,1))
707  qright(1) = half*(qrr(l,1)+qrl(l,1))
708  !
709  ! pressure
710  qleft(2) = half*(qll(l,2)+qlr(l,2))
711  qright(2) = half*(qrr(l,2)+qrl(l,2))
712  !
713  ! vt1 becomes normal velocity
714  qleft(3) = half*(qll(l,3)+qlr(l,3))
715  qright(3) = half*(qrr(l,3)+qrl(l,3))
716  !
717  ! bt1 becomes normal magnetic field
718  qleft(4) = half*(qll(l,6)+qlr(l,6))
719  qright(4) = half*(qrr(l,6)+qrl(l,6))
720  !
721  ! vt2 becomes transverse velocity field
722  qleft(5) = half*(qll(l,4)+qlr(l,4))
723  qright(5) = half*(qrr(l,4)+qrl(l,4))
724  !
725  ! bt2 becomes transverse magnetic field
726  qleft(6) = half*(qll(l,7)+qlr(l,7))
727  qright(6) = half*(qrr(l,7)+qrl(l,7))
728  !
729  ! velocity component perp. to the plane is now transverse
730  qleft(7) = half*(qll(l,5)+qlr(l,5))
731  qright(7) = half*(qrr(l,5)+qrl(l,5))
732  !
733  ! magnetic field component perp. to the plane is now transverse
734  qleft(8) = half*(qll(l,8)+qlr(l,8))
735  qright(8) = half*(qrr(l,8)+qrl(l,8))
736  !
737  CALL lax_friedrich(qleft,qright,fmean_x)
738  !
739  ! call the second solver in the y direction
740  ! density
741  qleft(1) = half*(qll(l,1)+qrl(l,1))
742  qright(1) = half*(qrr(l,1)+qlr(l,1))
743  !
744  ! pressure
745  qleft(2) = half*(qll(l,2)+qrl(l,2))
746  qright(2) = half*(qrr(l,2)+qlr(l,2))
747  !
748  ! vt2 becomes normal velocity
749  qleft(3) = half*(qll(l,4)+qrl(l,4))
750  qright(3) = half*(qrr(l,4)+qlr(l,4))
751  !
752  ! bt2 becomes normal magnetic field
753  qleft(4) = half*(qll(l,7)+qrl(l,7))
754  qright(4) = half*(qrr(l,7)+qlr(l,7))
755  !
756  ! vt1 becomes transverse velocity field
757  qleft(5) = half*(qll(l,3)+qrl(l,3))
758  qright(5) = half*(qrr(l,3)+qlr(l,3))
759  !
760  ! bt1 becomes transverse magnetic field
761  qleft(6) = half*(qll(l,6)+qrl(l,6))
762  qright(6) = half*(qrr(l,6)+qlr(l,6))
763  !
764  ! velocity component perp. to the plane is now transverse
765  qleft(7) = half*(qll(l,5)+qrl(l,5))
766  qright(7) = half*(qrr(l,5)+qlr(l,5))
767  !
768  ! magnetic field component perp. to the plane is now transverse
769  qleft(8) = half*(qll(l,8)+qrl(l,8))
770  qright(8) = half*(qrr(l,8)+qlr(l,8))
771  !
772  CALL lax_friedrich(qleft,qright,fmean_y)
773  !
774  e = ev + (fmean_x(6) - fmean_y(6))
775  if (xdim == 1) then
776  emf(l,j,k) = e
777  else if (xdim == 2) then
778  emf(k,l,j) = e
779  else if (xdim == 3) then
780  emf(j,k,l) = e
781  endif
782  end do
783  !call trace%end (itimer, detailed_timer=detailed_timer)
784 END SUBROUTINE llf_2d
785 
786 !===============================================================================
787 !> LLF
788 !===============================================================================
789 SUBROUTINE llf (qleft, qright, fgdnv, ngrid, detailed_timer)
790  real(DP), dimension(:,:):: qleft, qright, fgdnv
791  integer:: ngrid, l
792  logical, optional:: detailed_timer
793  !-----------------------------------------------------------------------------
794  do l=1,ngrid
795  call lax_friedrich (qleft(l,:), qright(l,:), fgdnv(l,:))
796  end do
797 END SUBROUTINE llf
798 
799 !===============================================================================
800 !> 1D local Lax-Friedrich Riemann solver
801 !===============================================================================
802 SUBROUTINE lax_friedrich (qleft, qright, fgdnv)
803  REAL(dp), DIMENSION(:):: qleft, qright, fgdnv
804  REAL(dp), DIMENSION(size(qleft)):: fleft, fright, fmean
805  REAL(dp), DIMENSION(size(qleft)):: uleft ,uright, udiff
806  REAL(dp):: vleft, vright, bx_mean
807  !-----------------------------------------------------------------------------
808  ! Enforce continuity of normal component
809  bx_mean = half*(qleft(4)+qright(4))
810  qleft(4) = bx_mean
811  qright(4) = bx_mean
812  !
813  CALL find_mhd_flux (qleft ,uleft ,fleft )
814  CALL find_mhd_flux (qright,uright,fright)
815  !
816  ! find the mean flux
817  fmean = half * ( fright + fleft )
818  !
819  ! find the largest eigenvalue in the normal direction to the interface
820  CALL find_speed_info (qleft ,vleft )
821  CALL find_speed_info (qright,vright)
822  !
823  ! difference between the 2 states
824  udiff = half * ( uright - uleft )
825  !
826  ! the local Lax-Friedrich flux
827  fgdnv = fmean - max(vleft,vright) * udiff
828 END SUBROUTINE lax_friedrich
829 
830 !===============================================================================
831 !===============================================================================
832 SUBROUTINE find_mhd_flux (qvar, cvar, ff)
833  !! compute the 1D MHD fluxes from the conservative variables
834  !! the structure of qvar is : rho, Pressure, Vnormal, Bnormal,
835  !! Vtransverse1, Btransverse1, Vtransverse2, Btransverse2
836  IMPLICIT NONE
837  INTEGER :: i , ivar, ib
838  REAL(dp), DIMENSION(:):: qvar, cvar, ff
839  REAL(dp):: ecin,emag,etot,d,u,v,w,a,b,c,p,ptot,entho
840  !-----------------------------------------------------------------------------
841  ! Local variables
842  entho = one/(gamma-one)
843  if (isothermal) then
844  entho = one
845  else
846  entho = one/max(gamma-one,1e-6)
847  end if
848  d=qvar(1); p=qvar(2); u=qvar(3); a=qvar(4)
849  v=qvar(5); b=qvar(6); w=qvar(7); c=qvar(8)
850  ecin = half*(u*u+v*v+w*w)*d
851  emag = half*(a*a+b*b+c*c)
852  etot = p*entho+ecin+emag
853  ptot = p + emag
854  !
855  ! Compute conservative variables
856  cvar(1) = d
857  cvar(2) = etot
858  cvar(3) = d*u
859  cvar(4) = a
860  cvar(5) = d*v
861  cvar(6) = b
862  cvar(7) = d*w
863  cvar(8) = c
864  !
865  ! Compute fluxes
866  ff(1) = d*u
867  ff(2) = (etot+ptot)*u-a*(a*u+b*v+c*w)
868  ff(3) = d*u*u+ptot-a*a
869  ff(4) = zero
870  ff(5) = d*u*v-a*b
871  ff(6) = b*u-a*v
872  ff(7) = d*u*w-a*c
873  ff(8) = c*u-a*w
874 END SUBROUTINE find_mhd_flux
875 
876 !===============================================================================
877 !===============================================================================
878 SUBROUTINE find_speed_fast(qvar,vel_info)
879  !! calculate the fast magnetosonic velocity
880  !! the structure of qvar is : rho, Pressure, Vnormal, Bnormal,
881  !! Vtransverse1,Btransverse1,Vtransverse2,Btransverse2
882  REAL(dp), DIMENSION(:):: qvar
883  REAL(dp) :: vel_info
884  REAL(dp) :: d,p,u,v,w,a,b,c,b2,c2,d2,cf
885  integer, save:: itimer=0
886  !-----------------------------------------------------------------------------
887  d=qvar(1); p=qvar(2); u=qvar(3); a=qvar(4)
888  v=qvar(5); b=qvar(6); w=qvar(7); c=qvar(8)
889  d=max(d,smallr)
890  p=max(p,smallr*smallc*smallc)
891  b2 = a*a+b*b+c*c
892  c2 = gamma*p/d
893  d2 = half*(b2/d+c2)
894  cf = sqrt( d2 + sqrt(max(d2*d2-c2*a*a/d,zero)) )
895  vel_info = cf
896 END SUBROUTINE find_speed_fast
897 
898 !===============================================================================
899 !===============================================================================
900 SUBROUTINE find_speed_info (qvar, vel_info)
901  !! calculate the fastest velocity at which information is exchanged
902  !! at the interface
903  !! the structure of qvar is : rho, Pressure, Vnormal, Bnormal,
904  !! Vtransverse1,Btransverse1,Vtransverse2,Btransverse2
905  IMPLICIT NONE
906  REAL(dp), DIMENSION(:):: qvar
907  REAL(dp):: vel_info
908  REAL(dp):: d, p, u, v, w, a, b, c, b2, c2, d2, cf
909  !-----------------------------------------------------------------------------
910  d=qvar(1); p=qvar(2); u=qvar(3); a=qvar(4)
911  v=qvar(5); b=qvar(6); w=qvar(7); c=qvar(8)
912  b2 = a*a+b*b+c*c
913  c2 = gamma*p/d
914  d2 = half*(b2/d+c2)
915  cf = sqrt( d2 + sqrt(d2**2-c2*a*a/d) )
916  vel_info = cf+abs(u)
917 END SUBROUTINE find_speed_info
918 
919 END MODULE riemann_mod
Definition: io_mod.f90:4