12 procedure,
nopass:: approx
13 procedure,
nopass:: acoustic
14 procedure,
nopass:: llf
15 procedure,
nopass:: hll
16 procedure,
nopass:: hllc
17 procedure,
nopass:: hllc_a
18 procedure,
nopass:: hllc_v
20 type(riemann_t),
public:: rieman
25 subroutine approx(qleft,qright,fgdnv,ngrid,nvar)
27 real(dp),
dimension(1:ngrid,1:nvar)::qleft,qright,qgdnv,fgdnv
30 real(dp),
dimension(1:ngrid)::rl ,ul ,pl ,cl
31 real(dp),
dimension(1:ngrid)::rr ,ur ,pr ,cr
32 real(dp),
dimension(1:ngrid)::ro ,uo ,po ,co
33 real(dp),
dimension(1:ngrid)::rstar,ustar,pstar,cstar
34 real(dp),
dimension(1:ngrid)::wl ,wr ,wo
35 real(dp),
dimension(1:ngrid)::sgnm ,spin ,spout,ushock
36 real(dp),
dimension(1:ngrid)::frac ,delp ,pold
37 integer ,
dimension(1:ngrid)::ind ,ind2
40 real(dp)::smallp, gamma6, ql, qr, usr, usl, wwl, wwr, smallpp, entho, etot
41 integer ::i, j, n, iter, n_new
44 smallp = smallc**2/gamma
45 smallpp = smallr*smallp
46 gamma6 = (gamma+one)/(two*gamma)
47 entho = one/(gamma-one)
51 rl(i)=max(qleft(i,1),smallr)
53 pl(i)=max(qleft(i,3),rl(i)*smallp)
54 rr(i)=max(qright(i,1),smallr)
56 pr(i)=max(qright(i,3),rr(i)*smallp)
61 cl(i) = gamma*pl(i)*rl(i)
62 cr(i) = gamma*pr(i)*rr(i)
67 wl(i)= sqrt(cl(i)); wr(i) = sqrt(cr(i))
68 pstar(i) = ((wr(i)*pl(i)+wl(i)*pr(i))+wl(i)*wr(i)*(ul(i)-ur(i)))/(wl(i)+wr(i))
69 pstar(i) = max(pstar(i),0.0_dp)
79 do iter = 1,niter_riemann
81 wwl=sqrt(cl(ind(i))*(one+gamma6*(pold(i)-pl(ind(i)))/pl(ind(i))))
82 wwr=sqrt(cr(ind(i))*(one+gamma6*(pold(i)-pr(ind(i)))/pr(ind(i))))
83 ql=two*wwl**3/(wwl**2+cl(ind(i)))
84 qr=two*wwr**3/(wwr**2+cr(ind(i)))
85 usl=ul(ind(i))-(pold(i)-pl(ind(i)))/wwl
86 usr=ur(ind(i))+(pold(i)-pr(ind(i)))/wwr
87 delp(i)=max(qr*ql/(qr+ql)*(usl-usr),-pold(i))
90 pold(i)=pold(i)+delp(i)
94 uo(i)=abs(delp(i)/(pold(i)+smallpp))
106 if(uo(i)<=1.d-06)
then 120 pstar(ind(i))=pold(i)
123 wl(i) = sqrt(cl(i)*(one+gamma6*(pstar(i)-pl(i))/pl(i)))
124 wr(i) = sqrt(cr(i)*(one+gamma6*(pstar(i)-pr(i))/pr(i)))
130 ustar(i) = half*(ul(i) + (pl(i)-pstar(i))/wl(i) + &
131 & ur(i) - (pr(i)-pstar(i))/wr(i) )
136 sgnm(i) = sign(one,ustar(i))
154 co(i) = max(smallc,sqrt(abs(gamma*po(i)/ro(i))))
159 if(pstar(i)>= po(i))
then 161 rstar(i) = ro(i)/(one+ro(i)*(po(i)-pstar(i))/wo(i)**2)
164 rstar(i) = ro(i)*(pstar(i)/po(i))**(one/gamma)
170 rstar(i) = max(rstar(i),smallr)
172 cstar(i) = sqrt(abs(gamma*pstar(i)/rstar(i)))
173 cstar(i) = max(cstar(i),smallc)
175 spout(i) = co(i)-sgnm(i)*uo(i)
176 spin(i) = cstar(i)-sgnm(i)*ustar(i)
178 ushock(i) = wo(i)/ro(i)-sgnm(i)*uo(i)
182 if(pstar(i)>=po(i))
then 190 if(spout(i)<=zero)
then 194 else if(spin(i)>=zero)
then 195 qgdnv(i,1) = rstar(i)
196 qgdnv(i,2) = ustar(i)
197 qgdnv(i,3) = pstar(i)
199 frac(i)=spout(i)/(spout(i)-spin(i))
200 qgdnv(i,2) = frac(i)*ustar(i) + (one - frac(i))*uo(i)
201 qgdnv(i,3) = frac(i)*pstar(i) + (one - frac(i))*po(i)
202 qgdnv(i,1) = ro(i)*(qgdnv(i,3)/po(i))**(one/gamma)
210 qgdnv(i,n) = qleft(i,n)
212 qgdnv(i,n) = qright(i,n)
219 fgdnv(i,1) = qgdnv(i,1)*qgdnv(i,2)
220 fgdnv(i,2) = qgdnv(i,3)+qgdnv(i,1)*qgdnv(i,2)**2
221 etot = qgdnv(i,3)*entho + half*qgdnv(i,1)*qgdnv(i,2)**2
222 etot = etot + half*qgdnv(i,1)*qgdnv(i,4)**2
223 etot = etot + half*qgdnv(i,1)*qgdnv(i,5)**2
224 fgdnv(i,3) = qgdnv(i,2)*(etot+qgdnv(i,3))
229 fgdnv(i,n) = fgdnv(i,1)*qgdnv(i,n)
232 end subroutine approx
236 subroutine acoustic(qleft,qright,fgdnv,ngrid,nvar)
238 real(dp),
dimension(1:ngrid,1:nvar)::qleft,qright,qgdnv,fgdnv
242 real(dp)::smallp, entho, etot
245 real(dp),
dimension(1:ngrid)::rl ,ul ,pl ,cl
246 real(dp),
dimension(1:ngrid)::rr ,ur ,pr ,cr
247 real(dp),
dimension(1:ngrid)::ro ,uo ,po ,co
248 real(dp),
dimension(1:ngrid)::rstar,ustar,pstar,cstar
249 real(dp),
dimension(1:ngrid)::wl ,wr ,wo
250 real(dp),
dimension(1:ngrid)::sgnm ,spin ,spout,ushock
251 real(dp),
dimension(1:ngrid)::frac
254 smallp = smallc**2/gamma
255 entho = one/(gamma-one)
259 rl(i)=max(qleft(i,1),smallr)
261 pl(i)=max(qleft(i,3),rl(i)*smallp)
262 rr(i)=max(qright(i,1),smallr)
264 pr(i)=max(qright(i,3),rr(i)*smallp)
269 cl(i) = sqrt(gamma*pl(i)/rl(i))
270 cr(i) = sqrt(gamma*pr(i)/rr(i))
273 pstar(i) = ((wr(i)*pl(i)+wl(i)*pr(i))+wl(i)*wr(i)*(ul(i)-ur(i))) &
275 ustar(i) = ((wr(i)*ur(i)+wl(i)*ul(i))+(pl(i)-pr(i))) &
282 sgnm(i) = sign(one,ustar(i))
304 rstar(i) = ro(i)+(pstar(i)-po(i))/co(i)**2
305 rstar(i) = max(rstar(i),smallr)
306 cstar(i) = sqrt(abs(gamma*pstar(i)/rstar(i)))
307 cstar(i) = max(cstar(i),smallc)
312 spout(i) = co(i)-sgnm(i)*uo(i)
313 spin(i) = cstar(i)-sgnm(i)*ustar(i)
318 ushock(i) = half*(spin(i)+spout(i))
319 ushock(i) = max(ushock(i),-sgnm(i)*ustar(i))
322 if(pstar(i)>=po(i))
then 330 if(spout(i)<zero)
then 334 else if(spin(i)>=zero)
then 335 qgdnv(i,1) = rstar(i)
336 qgdnv(i,2) = ustar(i)
337 qgdnv(i,3) = pstar(i)
339 frac(i) = spout(i)/(spout(i)-spin(i))
340 qgdnv(i,1) = frac(i)*rstar(i) + (one - frac(i))*ro(i)
341 qgdnv(i,2) = frac(i)*ustar(i) + (one - frac(i))*uo(i)
342 qgdnv(i,3) = frac(i)*pstar(i) + (one - frac(i))*po(i)
350 qgdnv(i,n) = qleft(i,n)
352 qgdnv(i,n) = qright(i,n)
359 fgdnv(i,1) = qgdnv(i,1)*qgdnv(i,2)
360 fgdnv(i,2) = qgdnv(i,3)+qgdnv(i,1)*qgdnv(i,2)**2
361 etot = qgdnv(i,3)*entho + half*qgdnv(i,1)*qgdnv(i,2)**2
362 etot = etot + half*qgdnv(i,1)*qgdnv(i,4)**2
363 etot = etot + half*qgdnv(i,1)*qgdnv(i,5)**2
364 fgdnv(i,3) = qgdnv(i,2)*(etot+qgdnv(i,3))
369 fgdnv(i,n) = fgdnv(i,1)*qgdnv(i,n)
372 end subroutine acoustic
376 subroutine llf(qleft,qright,fgdnv,ngrid,nvar)
378 real(dp),
dimension(1:ngrid,1:nvar)::qleft,qright,fgdnv
381 real(dp),
dimension(1:ngrid,1:nvar)::fleft,fright,uleft,uright
382 real(dp),
dimension(1:ngrid)::cmax
386 real(dp)::smallp, entho
387 real(dp)::rl ,ul ,pl ,cl
388 real(dp)::rr ,ur ,pr ,cr
391 smallp = smallc**2/gamma
392 entho = one/(gamma-one)
396 rl=max(qleft(i,1),smallr)
398 pl=max(qleft(i,3),rl*smallp)
399 rr=max(qright(i,1),smallr)
401 pr=max(qright(i,3),rr*smallp)
402 cl= sqrt(gamma*pl/rl)
403 cr= sqrt(gamma*pr/rr)
404 cmax(i)=max(abs(ul)+cl,abs(ur)+cr)
410 uleft(i,1) = qleft(i,1)
411 uright(i,1) = qright(i,1)
413 uleft(i,2) = qleft(i,1)*qleft(i,2)
414 uright(i,2) = qright(i,1)*qright(i,2)
416 uleft(i,3) = qleft(i,3)*entho + half*qleft(i,1)*qleft(i,2)**2
417 uright(i,3) = qright(i,3)*entho + half*qright(i,1)*qright(i,2)**2
418 uleft(i,3) = uleft(i,3) + half*qleft(i,1)*qleft(i,4)**2
419 uright(i,3) = uright(i,3) + half*qright(i,1)*qright(i,4)**2
420 uleft(i,3) = uleft(i,3) + half*qleft(i,1)*qleft(i,5)**2
421 uright(i,3) = uright(i,3) + half*qright(i,1)*qright(i,5)**2
426 uleft(i,n) = qleft(i,1)*qleft(i,n)
427 uright(i,n) = qright(i,1)*qright(i,n)
434 fleft(i,1) = uleft(i,2)
435 fright(i,1) = uright(i,2)
437 fleft(i,2) = qleft(i,3)+uleft(i,2)*qleft(i,2)
438 fright(i,2) = qright(i,3)+uright(i,2)*qright(i,2)
440 fleft(i,3) = qleft(i,2)*(uleft(i,3)+qleft(i,3))
441 fright(i,3) = qright(i,2)*(uright(i,3)+qright(i,3))
446 fleft(i,n) = fleft(i,1)*qleft(i,n)
447 fright(i,n) = fright(i,1)*qright(i,n)
454 fgdnv(i,n) = half*(fleft(i,n)+fright(i,n)-cmax(i)*(uright(i,n)-uleft(i,n)))
461 subroutine hllc (qleft,qright,fgdnv,ngrid,nvar,isothermal,error,detailed_timer)
463 real(dp),
dimension(1:ngrid,1:nvar)::qleft,qright,fgdnv
466 REAL(dp)::rl,pl,ul,ecinl,etotl,ptotl
467 REAL(dp)::rr,pr,ur,ecinr,etotr,ptotr
468 REAL(dp)::cfastl,rcl,rstarl
469 REAL(dp)::cfastr,rcr,rstarr
470 REAL(dp)::etotstarl,etotstarr
471 REAL(dp)::ustar,ptotstar
472 REAL(dp)::ro,uo,ptoto,etoto
475 LOGICAL :: isothermal,error,detailed_timer
476 INTEGER,
SAVE:: itimer=0
478 if (detailed_timer)
call trace_begin (
'riemann::hllc', itimer=itimer)
480 smallp = smallc**2/gamma
481 if (isothermal .or. gamma==1d0)
then 484 entho=one/(gamma-one)
488 if (io%verbose>1)
then 492 rl=max(qleft(i,1),smallr)
493 pl=max(qleft(i,3),rl*smallp)
496 ecinl = half*rl*(ul*ul+qleft(i,4)**2+qleft(i,5)**2)
497 etotl = pl*entho+ecinl
501 rr=max(qright(i,1),smallr)
502 pr=max(qright(i,3),rr*smallp)
505 ecinr = half*rr*(ur*ur+qright(i,4)**2+qright(i,5)**2)
506 etotr = pr*entho+ecinr
510 cfastl=sqrt(max(gamma*pl/rl,smallc**2))
511 cfastr=sqrt(max(gamma*pr/rr,smallc**2))
514 sl=min(ul,ur)-max(cfastl,cfastr)
515 sr=max(ul,ur)+max(cfastl,cfastr)
520 if ((rcl+rcl)==0.0)
then 521 print 1,
'HLLC: rl,Pl,ul,ecinl,etotl,Ptotl,cfastl,SL =', &
522 i,rl,pl,ul,ecinl,etotl,ptotl,cfastl,sl,sl-ul
523 print 1,
'HLLC: rr,Pr,ur,ecinr,etotr,Ptotr,cfastr,SR =', &
524 i,rr,pr,ur,ecinr,etotr,ptotr,cfastr,sr,sr-ur
525 1
format(a,i4,1p,10e12.3)
531 ustar =(rcr*ur +rcl*ul + (ptotl-ptotr))/(rcr+rcl)
532 ptotstar=(rcr*ptotl+rcl*ptotr+rcl*rcr*(ul-ur))/(rcr+rcl)
533 if ((sl-ustar)==0.0 .or. (sr-ustar)==0.0)
then 534 print 1,
'HLLC: rl,Pl,ul,ecinl,etotl,Ptotl,cfastl,SL =', &
535 i,rl,pl,ul,ecinl,etotl,ptotl,cfastl,sl,sl-ul
536 print 1,
'HLLC: rr,Pr,ur,ecinr,etotr,Ptotr,cfastr,SR =', &
537 i,rr,pr,ur,ecinr,etotr,ptotr,cfastr,sr,sr-ur
538 print 1,
'HLLC: ustar,Ptotstar', i,ustar,ptotstar
539 call mpi%abort (
'HLLC')
545 rstarl=rl*(sl-ul)/(sl-ustar)
546 etotstarl=((sl-ul)*etotl-ptotl*ul+ptotstar*ustar)/(sl-ustar)
549 rstarr=rr*(sr-ur)/(sr-ustar)
550 etotstarr=((sr-ur)*etotr-ptotr*ur+ptotstar*ustar)/(sr-ustar)
558 else if(ustar>0d0)
then 577 fgdnv(i,2) = ro*uo*uo+ptoto
578 fgdnv(i,3) = (etoto+ptoto)*uo
584 rl=max(qleft(i,1),smallr)
585 pl=max(qleft(i,3),rl*smallp)
588 ecinl = half*rl*(ul*ul+qleft(i,4)**2+qleft(i,5)**2)
589 etotl = pl*entho+ecinl
593 rr=max(qright(i,1),smallr)
594 pr=max(qright(i,3),rr*smallp)
597 ecinr = half*rr*(ur*ur+qright(i,4)**2+qright(i,5)**2)
598 etotr = pr*entho+ecinr
602 cfastl=sqrt(max(gamma*pl/rl,smallc**2))
603 cfastr=sqrt(max(gamma*pr/rr,smallc**2))
606 sl=min(ul,ur)-max(cfastl,cfastr)
607 sr=max(ul,ur)+max(cfastl,cfastr)
614 ustar =(rcr*ur +rcl*ul + (ptotl-ptotr))/(rcr+rcl)
615 ptotstar=(rcr*ptotl+rcl*ptotr+rcl*rcr*(ul-ur))/(rcr+rcl)
618 rstarl=rl*(sl-ul)/(sl-ustar)
619 etotstarl=((sl-ul)*etotl-ptotl*ul+ptotstar*ustar)/(sl-ustar)
622 rstarr=rr*(sr-ur)/(sr-ustar)
623 etotstarr=((sr-ur)*etotr-ptotr*ur+ptotstar*ustar)/(sr-ustar)
632 ro =merge(rstarr ,ro ,sr>0d0)
633 uo =merge(ustar ,uo ,sr>0d0)
634 ptoto=merge(ptotstar ,ptoto ,sr>0d0)
635 etoto=merge(etotstarr,etoto ,sr>0d0)
636 ro =merge(rstarr ,ro ,ustar>0d0)
637 uo =merge(ustar ,uo ,ustar>0d0)
638 ptoto=merge(ptotstar ,ptoto ,ustar>0d0)
639 etoto=merge(etotstarl,etoto ,ustar>0d0)
640 ro =merge(rl ,ro ,sl>0d0)
641 uo =merge(ul ,uo ,sl>0d0)
642 ptoto=merge(ptotl ,ptoto ,sl>0d0)
643 etoto=merge(etotl ,etoto ,sl>0d0)
650 else if(ustar>0d0)
then 670 fgdnv(i,2) = ro*uo*uo+ptoto
671 fgdnv(i,3) = (etoto+ptoto)*uo
677 fgdnv(i,ivar) = fgdnv(i,1)*merge(qleft(i,ivar),qright(i,ivar),fgdnv(i,1)>0)
680 if (detailed_timer)
call trace_end (itimer)
685 subroutine hllc_a (qleft,qright,fgdnv,ngrid,nvar)
687 real(dp),
dimension(1:ngrid,1:nvar)::qleft,qright,fgdnv
694 REAL(dp),
dimension(ngrid)::sl,sr
695 REAL(dp),
dimension(ngrid)::ro,uo,ptoto,etoto
696 REAL(dp),
dimension(ngrid)::rl,ul,ptotl,etotl
697 REAL(dp),
dimension(ngrid)::rr,ur,ptotr,etotr
698 REAL(dp),
dimension(ngrid)::rstarl,etotstarl
699 REAL(dp),
dimension(ngrid)::rstarr,etotstarr
700 REAL(dp),
dimension(ngrid)::ustar,ptotstar
704 smallp = smallc**2/gamma
705 entho = one/(gamma-one)
710 rl(i)=max(qleft(i,1),smallr)
711 pl=max(qleft(i,3),rl(i)*smallp)
714 ecinl = half*rl(i)*(ul(i)*ul(i)+qleft(i,4)**2+qleft(i,5)**2)
715 etotl(i) = pl*entho+ecinl
719 rr(i)=max(qright(i,1),smallr)
720 pr=max(qright(i,3),rr(i)*smallp)
723 ecinr = half*rr(i)*(ur(i)*ur(i)+qright(i,4)**2+qright(i,5)**2)
724 etotr(i) = pr*entho+ecinr
728 cfastl=sqrt(max(gamma*pl/rl(i),smallc**2))
729 cfastr=sqrt(max(gamma*pr/rr(i),smallc**2))
732 sl(i)=min(ul(i),ur(i))-max(cfastl,cfastr)
733 sr(i)=max(ul(i),ur(i))+max(cfastl,cfastr)
736 rcl=rl(i)*(ul(i)-sl(i))
737 rcr=rr(i)*(sr(i)-ur(i))
740 ustar(i) =(rcr*ur(i) +rcl*ul(i) + (ptotl(i)-ptotr(i)))/(rcr+rcl)
741 ptotstar(i)=(rcr*ptotl(i)+rcl*ptotr(i)+rcl*rcr*(ul(i)-ur(i)))/(rcr+rcl)
746 rstarl(i)=rl(i)*(sl(i)-ul(i))/(sl(i)-ustar(i))
747 etotstarl(i)=((sl(i)-ul(i))*etotl(i)-ptotl(i)*ul(i)+ptotstar(i)*ustar(i))/(sl(i)-ustar(i))
750 rstarr(i)=rr(i)*(sr(i)-ur(i))/(sr(i)-ustar(i))
751 etotstarr(i)=((sr(i)-ur(i))*etotr(i)-ptotr(i)*ur(i)+ptotstar(i)*ustar(i))/(sr(i)-ustar(i))
761 else if(ustar(i)>0d0)
then 765 etoto(i)=etotstarl(i)
766 else if (sr(i)>0d0)
then 770 etoto(i)=etotstarr(i)
781 fgdnv(i,1) = ro(i)*uo(i)
782 fgdnv(i,2) = ro(i)*uo(i)*uo(i)+ptoto(i)
783 fgdnv(i,3) = (etoto(i)+ptoto(i))*uo(i)
788 fgdnv(i,ivar) = merge(fgdnv(i,1)*qleft(i,ivar), &
789 fgdnv(i,1)*qright(i,ivar), fgdnv(i,1)>0)
792 end subroutine hllc_a
796 subroutine hllc_v (qleft,qright,fgdnv,ngrid,nvar)
798 real(dp),
dimension(1:ngrid,1:nvar)::qleft,qright,fgdnv
799 REAL(dp),
dimension(ngrid)::sl,sr
800 REAL(dp),
dimension(ngrid)::entho
801 REAL(dp),
dimension(ngrid)::rl,pl,ul,ecinl,etotl,ptotl
802 REAL(dp),
dimension(ngrid)::rr,pr,ur,ecinr,etotr,ptotr
803 REAL(dp),
dimension(ngrid)::cfastl,rcl,rstarl
804 REAL(dp),
dimension(ngrid)::cfastr,rcr,rstarr
805 REAL(dp),
dimension(ngrid)::etotstarl,etotstarr
806 REAL(dp),
dimension(ngrid)::ustar,ptotstar
807 REAL(dp),
dimension(ngrid)::ro,uo,ptoto,etoto
808 REAL(dp),
dimension(ngrid)::smallp
812 smallp = smallc**2/gamma
813 entho = one/(gamma-one)
817 rl=max(qleft(:,1),smallr)
818 pl=max(qleft(:,3),rl*smallp)
821 ecinl = half*rl*(ul*ul+qleft(:,4)**2+qleft(:,5)**2)
822 etotl = pl*entho+ecinl
826 rr=max(qright(:,1),smallr)
827 pr=max(qright(:,3),rr*smallp)
830 ecinr = half*rr*(ur*ur+qright(:,4)**2+qright(:,5)**2)
831 etotr = pr*entho+ecinr
835 cfastl=sqrt(max(gamma*pl/rl,smallc**2))
836 cfastr=sqrt(max(gamma*pr/rr,smallc**2))
839 sl=min(ul,ur)-max(cfastl,cfastr)
840 sr=max(ul,ur)+max(cfastl,cfastr)
847 ustar =(rcr*ur +rcl*ul + (ptotl-ptotr))/(rcr+rcl)
848 ptotstar=(rcr*ptotl+rcl*ptotr+rcl*rcr*(ul-ur))/(rcr+rcl)
851 rstarl=rl*(sl-ul)/(sl-ustar)
852 etotstarl=((sl-ul)*etotl-ptotl*ul+ptotstar*ustar)/(sl-ustar)
855 rstarr=rr*(sr-ur)/(sr-ustar)
856 etotstarr=((sr-ur)*etotr-ptotr*ur+ptotstar*ustar)/(sr-ustar)
864 else where (ustar>0d0)
883 fgdnv(:,2) = ro*uo*uo+ptoto
884 fgdnv(:,3) = (etoto+ptoto)*uo
887 fgdnv(:,ivar) = fgdnv(:,1)*qleft(:,ivar)
889 fgdnv(:,ivar) = fgdnv(:,1)*qright(:,ivar)
892 end subroutine hllc_v
896 subroutine hll(qleft,qright,fgdnv,ngrid,nvar)
898 real(dp),
dimension(1:ngrid,1:nvar)::qleft,qright,fgdnv
899 real(dp),
dimension(1:ngrid,1:nvar)::fleft,fright,uleft,uright
900 real(dp),
dimension(1:ngrid)::sl,sr
902 real(dp)::smallp, entho
903 real(dp)::rl ,ul ,pl ,cl
904 real(dp)::rr ,ur ,pr ,cr
907 smallp = smallc**2/gamma
908 entho = one/(gamma-one)
912 rl=max(qleft(i,1),smallr)
914 pl=max(qleft(i,3),rl*smallp)
915 rr=max(qright(i,1),smallr)
917 pr=max(qright(i,3),rr*smallp)
918 cl= sqrt(gamma*pl/rl)
919 cr= sqrt(gamma*pr/rr)
920 sl(i)=min(min(ul,ur)-max(cl,cr),zero)
921 sr(i)=max(max(ul,ur)+max(cl,cr),zero)
926 uleft(i,1) = qleft(i,1)
927 uright(i,1) = qright(i,1)
928 uleft(i,2) = qleft(i,1)*qleft(i,2)
929 uright(i,2) = qright(i,1)*qright(i,2)
930 uleft(i,3) = qleft(i,3)*entho + half*qleft(i,1)*qleft(i,2)**2
931 uright(i,3) = qright(i,3)*entho + half*qright(i,1)*qright(i,2)**2
932 uleft(i,3) = uleft(i,3) + half*qleft(i,1)*qleft(i,4)**2
933 uright(i,3) = uright(i,3) + half*qright(i,1)*qright(i,4)**2
934 uleft(i,3) = uleft(i,3) + half*qleft(i,1)*qleft(i,5)**2
935 uright(i,3) = uright(i,3) + half*qright(i,1)*qright(i,5)**2
940 uleft(i,n) = qleft(i,1)*qleft(i,n)
941 uright(i,n) = qright(i,1)*qright(i,n)
947 fleft(i,1) = uleft(i,2)
948 fright(i,1) = uright(i,2)
949 fleft(i,2) = qleft(i,3)+uleft(i,2)*qleft(i,2)
950 fright(i,2) = qright(i,3)+uright(i,2)*qright(i,2)
951 fleft(i,3) = qleft(i,2)*(uleft(i,3)+qleft(i,3))
952 fright(i,3) = qright(i,2)*(uright(i,3)+qright(i,3))
957 fleft(i,n) = fleft(i,1)*qleft(i,n)
958 fright(i,n) = fright(i,1)*qright(i,n)
965 fgdnv(i,n) = (sr(i)*fleft(i,n)-sl(i)*fright(i,n) &
966 & + sr(i)*sl(i)*(uright(i,n)-uleft(i,n)))/(sr(i)-sl(i))
971 end module riemann_mod