subroutine HLLC(leftst,rigtst,nflux)
!$acc routine seq
!=====================================================================
!
! HLLC Scheme
!
! Purpose
! Calculation of Numerical Flux by HLLC method
!
! Reference
! Toro EF, Spruce M, Speares W. (1992,1994)
!
! Input
! Output
!=====================================================================
use fluxmod, only: mflx,madd &
& , mudn,muvu,muvv,muvw,muet &
& , mfdn,mfvu,mfvv,mfvw,mfet &
& , mcsp,mvel,mpre &
& , mden,mrvu,mrvv,mrvw,meto &
& , mubu,mubv,mubw,mubp &
& , mfbu,mfbv,mfbw,mfbp &
& , mbmu,mbmv,mbmw &
& , ncomp, mst, med
implicit none
real(8),dimension(2*mflx+madd),intent(in)::leftst,rigtst
real(8),dimension(mflx),intent(out)::nflux
!----- U -----
! qql :: left state
! qqr :: right state
real(8) :: rol,vxl,vyl,vzl,ptl,eel
real(8) :: ror,vxr,vyr,vzr,ptr,eer
real(8) :: rxl,ryl,rzl
real(8) :: rxr,ryr,rzr
real(8) :: byl,bzl
real(8) :: byr,bzr
real(8) :: ptst
real(8),dimension(ncomp) :: scl,scr
!----- U* ----
! qqlst :: left state
! qqrst :: right state
real(8) :: rolst,vxlst,vylst,vzlst,eelst
real(8) :: rorst,vxrst,vyrst,vzrst,eerst
real(8) :: rxlst,rylst,rzlst
real(8) :: rxrst,ryrst,rzrst
real(8) :: bylst,bzlst
real(8) :: byrst,bzrst
real(8),dimension(ncomp) :: sclst,scrst
!----- flux ---
! fqql :: left physical flux
! fqqr :: right physical flux
real(8) :: frol,frxl,fryl,frzl,feel
real(8) :: fbyl,fbzl
real(8) :: fror,frxr,fryr,frzr,feer
real(8) :: fbyr,fbzr
real(8),dimension(ncomp) :: fscl,fscr
!----- wave speed ---
! sl :: left-going fastest signal velocity
! sr :: right-going fastest signal velocity
! sm :: contact discontinuity velocity
! slst :: left-going alfven velocity
! srst :: right-going alfven velocity
real(8) :: sm,sl,sr
! cfl :: left-state Fast wave velocity
! cfr :: right-sate Fast wave velocity
real(8) :: cfl,cfr
!--------------------
! temporary variables
real(8) :: sdl,sdr,sdml,sdmr,isdml,isdmr,rosdl,rosdr
real(8) :: temp
real(8),dimension(ncomp) :: scsdl,scsdr
! no if
real(8) :: sign1,maxs1,mins1
real(8) :: msl,msr
!----- Step 0. ----------------------------------------------------------|
!---- Left state
rol = leftst(mudn)
eel = leftst(muet)
rxl = leftst(muvu)
ryl = leftst(muvv)
rzl = leftst(muvw)
vxl = leftst(muvu)/leftst(mudn)
vyl = leftst(muvv)/leftst(mudn)
vzl = leftst(muvw)/leftst(mudn)
scl(1:ncomp) = leftst(mst:med)
byl = leftst(mubv)
bzl = leftst(mubw)
ptl = leftst(mpre)
!---- Right state
ror = rigtst(mudn)
eer = rigtst(muet)
rxr = rigtst(muvu)
ryr = rigtst(muvv)
rzr = rigtst(muvw)
vxr = rigtst(muvu)/rigtst(mudn)
vyr = rigtst(muvv)/rigtst(mudn)
vzr = rigtst(muvw)/rigtst(mudn)
scr(1:ncomp) = rigtst(mst:med)
byr = rigtst(mubv)
bzr = rigtst(mubw)
ptr = rigtst(mpre)
!----- Step 1. ----------------------------------------------------------|
! Compute wave left & right wave speed
!
cfl = leftst(mcsp)
cfr = rigtst(mcsp)
sl = min(vxl,vxr)-max(cfl,cfr) ! note sl is negative
sr = max(vxl,vxr)+max(cfl,cfr)
!----- Step 2. ----------------------------------------------------------|
! compute L/R fluxs
!
! Left value
frol = leftst(mfdn)
feel = leftst(mfet)
frxl = leftst(mfvu)
fryl = leftst(mfvv)
frzl = leftst(mfvw)
fbyl = leftst(mfbv)
fbzl = leftst(mfbw)
fscl(1:ncomp) = leftst(mflx+mst:mflx+med)
! Right value
! Left value
fror = rigtst(mfdn)
feer = rigtst(mfet)
frxr = rigtst(mfvu)
fryr = rigtst(mfvv)
frzr = rigtst(mfvw)
fbyr = rigtst(mfbv)
fbzr = rigtst(mfbw)
fscr(1:ncomp) = rigtst(mflx+mst:mflx+med)
!----- Step 4. ----------------------------------------------------------|
! compute middle and alfven wave
!
sdl = sl - vxl
sdr = sr - vxr
rosdl = rol*sdl
rosdr = ror*sdr
scsdl(1:ncomp) = scl(1:ncomp)*sdl
scsdr(1:ncomp) = scr(1:ncomp)*sdr
temp = 1.0d0/(rosdr - rosdl)
! Eq. 45
sm = (rosdr*vxr - rosdl*vxl - ptr + ptl)*temp
sdml = sl - sm; isdml = 1.0d0/sdml
sdmr = sr - sm; isdmr = 1.0d0/sdmr
!----- Step 5. ----------------------------------------------------------|
! compute intermediate states
!
! Eq. 49
ptst = (rosdr*ptl-rosdl*ptr+rosdl*rosdr*(vxr-vxl))*temp
!----- Step 5A. ----------------------------------------------------------|
! compute Ul*
!
rolst = rol*sdl *isdml
vxlst = sm
rxlst = rolst*vxlst
vylst = vyl
rylst = rolst*vylst
vzlst = vzl
rzlst = rolst*vzlst
bylst = rolst/rol * byl
bzlst = rolst/rol * bzl
eelst =(sdl*eel - ptl*vxl + ptst*sm )*isdml
!----- Step 5B. ----------------------------------------------------------|
! compute Ur*
!
rorst = rosdr *isdmr
vxrst = sm
rxrst = rorst*vxrst
vyrst = vyr
ryrst = rorst*vyrst
vzrst = vzr
rzrst = rorst*vzrst
byrst = rorst/ror * byr
bzrst = rorst/ror * bzr
eerst = (sdr*eer - ptr*vxr + ptst*sm )*isdmr
!----- Step 6. ----------------------------------------------------------|
! compute flux
sign1 = sign(1.0d0,sm) ! 1 for sm>0, -1 for sm<0
maxs1 = max(0.0d0,sign1) ! 1 sm>0, 0 for sm<0
mins1 = -min(0.0d0,sign1) ! 0 sm>0,-1 for sm<0
msl = min(sl ,0.0d0) ! 0 for sl > 0, sl for sl < 0
msr = max(sr ,0.0d0) ! S_R > 0
nflux(mden) = (frol+msl*(rolst-rol))*maxs1 &
& +(fror+msr*(rorst-ror))*mins1
nflux(meto) = (feel+msl*(eelst-eel))*maxs1 &
& +(feer+msr*(eerst-eer))*mins1
nflux(mrvu) = (frxl+msl*(rxlst-rxl))*maxs1 &
& +(frxr+msr*(rxrst-rxr))*mins1
nflux(mrvv) = (fryl+msl*(rylst-ryl))*maxs1 &
& +(fryr+msr*(ryrst-ryr))*mins1
nflux(mrvw) = (frzl+msl*(rzlst-rzl))*maxs1 &
& +(frzr+msr*(rzrst-rzr))*mins1
nflux(mbmu) = 0.0d0
nflux(mbmv) = (fbyl+msl*(bylst-byl))*maxs1 &
& +(fbyr+msr*(byrst-byr))*mins1
nflux(mbmw) = (fbzl+msl*(bzlst-bzl))*maxs1 &
& +(fbzr+msr*(bzrst-bzr))*mins1
nflux(mst:med) = &
& (fscl(1:ncomp)+msl*(sclst(1:ncomp)-scl(1:ncomp)))*maxs1 &
& +(fscr(1:ncomp)+msr*(scrst(1:ncomp)-scr(1:ncomp)))*mins1
return
end subroutine HLLC