mod_mhd_roe.t

Go to the documentation of this file.

!> Subroutines for Roe-type Riemann solver for MHD
module mod_mhd_roe
  use mod_mhd_phys
  use mod_functions_bfield, only: mag
  use mod_physics_roe
 
  implicit none
  private
 
  integer, parameter :: fastRW_ = 3,fastlw_=4,slowrw_=5,slowlw_=6 ! Characteristic
  integer, parameter :: entroW_ = 8,diverw_=7,alfvrw_=1,alfvlw_=2 ! waves
 
  public :: mhd_roe_init
 
contains
 
  subroutine mhd_roe_init()
    use mod_global_parameters, only: entropycoef, nw
    integer :: il
 
    phys_average => mhd_average
    phys_get_eigenjump => mhd_get_eigenjump
    phys_rtimes => mhd_rtimes
 
    nworkroe=15
    allocate(entropycoef(nw))
 
    do il = 1, nw
       select case(il)
       case(fastrw_,fastlw_,slowrw_,slowlw_)
          entropycoef(il) = 0.2d0
       case(alfvrw_,alfvlw_)
          entropycoef(il) = 0.4d0
       case default
          entropycoef(il) = -1.0d0
       end select
    end do
 
  end subroutine mhd_roe_init
 
  ! Eight-wave MHD Riemann solver. See Powell, Notes on the eigensystem, Gombosi
  ! Calculate the wroe average of primitive variables in wL and wR, assignment:
  ! rho -> sqrho, m -> v, e -> p, B_idim -> B_idim, B_idir -> beta_idir
  ! Calculate also alpha_f,alpha_s,c_f,c_s,csound2,dp,rhodv
  !
  ! wL,wR,wroe are all interface centered quantities
  subroutine mhd_average(wL,wR,x,ix^L,idim,wroe,workroe)
    use mod_global_parameters
 
    integer, intent(in)             :: ix^L, idim
    double precision, intent(in)    :: wL(ixG^T, nw), wR(ixG^T, nw)
    double precision, intent(inout) :: wroe(ixG^T, nw)
    double precision, intent(inout) :: workroe(ixG^T, nworkroe)
    double precision, intent(in)    :: x(ixG^T, 1:^ND)
 
    call average2(wl,wr,x,ix^l,idim,wroe,workroe(ixg^t,1),workroe(ixg^t,2), &
         workroe(ixg^t,3),workroe(ixg^t,4),workroe(ixg^t,5),workroe(ixg^t,6), &
         workroe(ixg^t,7),workroe(ixg^t,8))
 
  end subroutine mhd_average
 
  ! Eight-wave MHD Riemann solver. See Powell, Notes on the eigensystem, Gombosi
  ! Calculate the wroe average of primitive variables in wL and wR, assignment:
  ! rho -> sqrho, m -> v, e -> p, B_idim -> B_idim, B_idir -> beta_idir
  ! Calculate also alpha_f,alpha_s,c_f,c_s,csound2,dp,rhodv
  !
  ! wL,wR,wroe are all interface centered quantities
  subroutine average2(wL,wR,x,ix^L,idim,wroe,cfast,cslow,afast,aslow,csound2,dp, &
       rhodv,tmp)
    use mod_global_parameters
 
    integer                               :: ix^L,idim,idir,jdir,iw
    double precision, dimension(ixG^T,nw) :: wL,wR,wroe
    double precision, intent(in)          :: x(ixG^T,1:ndim)
    double precision, dimension(ixG^T)    :: cfast,cslow,afast,aslow,csound2,dp, &
         rhodv,tmp
 
    if (ndir==1) call mpistop("MHD with d=11 is the same as HD")
 
    !Averaging primitive variables
    wroe(ix^s,rho_)=half*(wl(ix^s,rho_)+wr(ix^s,rho_))
    do idir=1,ndir
      wroe(ix^s,mom(idir))=half*(wl(ix^s,mom(idir))/wl(ix^s,rho_)+wr(ix^s,mom(idir))/wr(ix^s,rho_))
      wroe(ix^s,mag(idir))=half*(wl(ix^s,mag(idir))+wr(ix^s,mag(idir)))
    end do
    ! Use afast and aslow for pressures pL and pR
    call mhd_get_pthermal(wl,x,ixg^ll,ix^l,afast)
    call mhd_get_pthermal(wr,x,ixg^ll,ix^l,aslow)
 
    if(mhd_energy) then
      wroe(ix^s,e_)=half*(afast(ix^s)+aslow(ix^s))
      ! dp=pR-pL
      dp(ix^s)=aslow(ix^s)-afast(ix^s)
    else
      dp(ix^s)=aslow(ix^s)-afast(ix^s)
    end if
 
    !CONSERVATIVE rho*dv_idim=dm_idim-v_idim*drho
    rhodv(ix^s)=wr(ix^s,mom(idim))-wl(ix^s,mom(idim))-&
         wroe(ix^s,mom(idim))*(wr(ix^s,rho_)-wl(ix^s,rho_))
 
    !Calculate csound2,cfast,cslow,alphafast and alphaslow
 
    ! get csound**2
    if(mhd_energy) then
      csound2(ix^s)=mhd_gamma*wroe(ix^s,p_)/wroe(ix^s,rho_)
    else
      csound2(ix^s)=mhd_gamma*mhd_adiab*wroe(ix^s,rho_)**(mhd_gamma-one)
    end if
 
    ! aa=B**2/rho+a**2
    cfast(ix^s)=sum(wroe(ix^s,mag(:))**2,dim=ndim+1)/wroe(ix^s,rho_)+csound2(ix^s)
 
    ! cs**2=0.5*(aa+dsqrt(aa**2-4*a**2*(b_i**2/rho)))
    cslow(ix^s)=half*(cfast(ix^s)-dsqrt(cfast(ix^s)**2-&
         4d0*csound2(ix^s)*wroe(ix^s,mag(idim))**2/wroe(ix^s,rho_)))
 
    ! cf**2=aa-cs**2
    cfast(ix^s)=cfast(ix^s)-cslow(ix^s)
 
    ! alpha_f**2=(a**2-cs**2)/(cf**2-cs**2)
    afast(ix^s)=(csound2(ix^s)-cslow(ix^s))/(cfast(ix^s)-cslow(ix^s))
    afast(ix^s)=min(one,max(afast(ix^s),zero))
 
    ! alpha_s=dsqrt(1-alpha_f**2)
    aslow(ix^s)=dsqrt(one-afast(ix^s))
 
    ! alpha_f=dsqrt(alpha_f**2)
    afast(ix^s)=dsqrt(afast(ix^s))
 
    ! cf=dsqrt(cf**2)
    cfast(ix^s)=dsqrt(cfast(ix^s))
 
    ! cs=dsqrt(cs**2)
    cslow(ix^s)=dsqrt(cslow(ix^s))
 
    !Replace the primitive variables with more useful quantities:
    ! rho -> dsqrt(rho)
    wroe(ix^s,rho_)=dsqrt(wroe(ix^s,rho_))
 
    ! Avoid sgn(b_idim)==0
    where(dabs(wroe(ix^s,mag(idim)))<smalldouble)&
         wroe(ix^s,mag(idim))=smalldouble
    ! B_idir,jdir -> beta_idir,jdir
    idir=idim+1-ndir*(idim/ndir)
    if(ndir==2)then
       where(wroe(ix^s,mag(idir))>=zero)
          wroe(ix^s,mag(idir))=one
       elsewhere
          wroe(ix^s,mag(idir))=-one
       end where
    else
       !beta_j=B_j/dsqrt(B_i**2+B_j**2); beta_i=B_i/dsqrt(B_i**2+B_j**2)
       jdir=idir+1-ndir*(idir/ndir)
       tmp(ix^s)=dsqrt(wroe(ix^s,mag(idir))**2+wroe(ix^s,mag(jdir))**2)
       where(tmp(ix^s)>smalldouble)
          wroe(ix^s,mag(idir))=wroe(ix^s,mag(idir))/tmp(ix^s)
          wroe(ix^s,mag(jdir))=wroe(ix^s,mag(jdir))/tmp(ix^s)
       elsewhere
          wroe(ix^s,mag(idir))=dsqrt(half)
          wroe(ix^s,mag(jdir))=dsqrt(half)
       end where
    endif
 
  end subroutine average2
 
  ! Calculate the il-th characteristic speed and the jump in the il-th
  ! characteristic variable in the idim direction within ixL.
  ! The eigenvalues and the l=r**(-1) matrix is calculated from wroe.
  ! jump(il)=Sum_il l(il,iw)*(wR(iw)-wL(iw)), where w are the conservative
  ! variables. However part of the summation is done in advance and saved into
  ! bdv,bdb,dp and dv variables. "smalla" contains a lower limit for "a" to be
  ! used in the entropy fix.
  !
  ! All the variables are centered on the cell interface, thus the 
  ! "*C" notation is omitted for sake of brevity.
  subroutine mhd_get_eigenjump(wL,wR,wroe,x,ix^L,il,idim,smalla,a,jump,workroe)
    use mod_global_parameters
 
    integer, intent(in)                         :: ix^L,il,idim
    double precision, dimension(ixG^T,nw)       :: wL,wR,wroe
    double precision, intent(in)                :: x(ixG^T,1:ndim)
    double precision, dimension(ixG^T)          :: smalla,a,jump
    double precision, dimension(ixG^T,nworkroe) :: workroe
 
    call geteigenjump2(wl,wr,wroe,x,ix^l,il,idim,smalla,a,jump, &
         workroe(ixg^t,1),workroe(ixg^t,2), &
         workroe(ixg^t,3),workroe(ixg^t,4),workroe(ixg^t,5),workroe(ixg^t,6), &
         workroe(ixg^t,7),workroe(ixg^t,8),workroe(ixg^t,9),workroe(ixg^t,10), &
         workroe(ixg^t,11),workroe(ixg^t,12),workroe(ixg^t,13))
 
  end subroutine mhd_get_eigenjump
 
  ! Calculate the il-th characteristic speed and the jump in the il-th
  ! characteristic variable in the idim direction within ixL.
  ! The eigenvalues and the l=r**(-1) matrix is calculated from wroe.
  ! jump(il)=Sum_il l(il,iw)*(wR(iw)-wL(iw)), where w are the conservative
  ! variables. However part of the summation is done in advance and saved into
  ! bdv,bdb,dp and dv variables. "smalla" contains a lower limit for "a" to be
  ! used in the entropy fix.
  !
  ! All the variables are centered on the cell interface, thus the 
  ! "*C" notation is omitted for sake of brevity.
  subroutine geteigenjump2(wL,wR,wroe,x,ix^L,il,idim,smalla,a,jump, &
       cfast,cslow,afast,aslow,csound2,dp,rhodv,bdv,bdb,cs2L,cs2R,cs2ca2L,cs2ca2R)
    use mod_global_parameters
    use mod_tvd
 
    integer                               :: ix^L,il,idim,idir,jdir
    double precision, dimension(ixG^T,nw) :: wL,wR,wroe
    double precision, intent(in)          :: x(ixG^T,1:ndim)
    double precision, dimension(ixG^T)    :: smalla,a,jump
    double precision, dimension(ixG^T)    :: cfast,cslow,afast,aslow,csound2,dp,rhodv
    double precision, dimension(ixG^T)    :: bdv,bdb
    double precision, dimension(ixG^T)    :: aL,aR,cs2L,cs2R,cs2ca2L,cs2ca2R
 
    idir=idim+1-ndir*(idim/ndir)
    jdir=idir+1-ndir*(idir/ndir)
 
    if(il==fastrw_)then
       !Fast and slow waves use bdv=sqrho**2*sign(bx)*(betay*dvy+betaz*dvz)
       !                        bdb=sqrho*a*          (betay*dBy+betaz*dBz)
       bdv(ix^s)=wroe(ix^s,mag(idir))* &
            (wr(ix^s,mom(idir))/wr(ix^s,rho_)-wl(ix^s,mom(idir))/wl(ix^s,rho_))
       if(ndir==3)bdv(ix^s)=bdv(ix^s)+wroe(ix^s,mag(jdir))* &
            (wr(ix^s,mom(jdir))/wr(ix^s,rho_)-wl(ix^s,mom(jdir))/wl(ix^s,rho_))
       bdv(ix^s)=bdv(ix^s)*sign(wroe(ix^s,rho_)**2,wroe(ix^s,mag(idim)))
 
       bdb(ix^s)=wroe(ix^s,mag(idir))*(wr(ix^s,mag(idir))-wl(ix^s,mag(idir)))
       if(ndir==3)bdb(ix^s)=bdb(ix^s)+&
            wroe(ix^s,mag(jdir))*(wr(ix^s,mag(jdir))-wl(ix^s,mag(jdir)))
       bdb(ix^s)=bdb(ix^s)*dsqrt(csound2(ix^s))*wroe(ix^s,rho_)
    endif
 
    if(il==alfvrw_)then
       !Alfven waves use      bdv=0.5*sqrho**2*      (betaz*dvy-betay*dvz)
       !                      bdb=0.5*sqrho*sign(bx)*(betaz*dBy-betay*dBz)
       bdv(ix^s)=wroe(ix^s,mag(jdir))* &
            (wr(ix^s,mom(idir))/wr(ix^s,rho_)-wl(ix^s,mom(idir))/wl(ix^s,rho_)) &
            -wroe(ix^s,mag(idir))* &
            (wr(ix^s,mom(jdir))/wr(ix^s,rho_)-wl(ix^s,mom(jdir))/wl(ix^s,rho_))
       bdb(ix^s)=wroe(ix^s,mag(jdir))*(wr(ix^s,mag(idir))-wl(ix^s,mag(idir))) &
            -wroe(ix^s,mag(idir))*(wr(ix^s,mag(jdir))-wl(ix^s,mag(jdir)))
       bdv(ix^s)=bdv(ix^s)*half*wroe(ix^s,rho_)**2
       bdb(ix^s)=bdb(ix^s)*half*sign(wroe(ix^s,rho_),wroe(ix^s,mag(idim)))
    endif
 
    select case(il)
    case(fastrw_)
       a(ix^s)=wroe(ix^s,mom(idim))+cfast(ix^s)
       jump(ix^s)=half/csound2(ix^s)*(&
            afast(ix^s)*(+cfast(ix^s)*rhodv(ix^s)+dp(ix^s))&
            +aslow(ix^s)*(-cslow(ix^s)*bdv(ix^s)+bdb(ix^s)))
    case(fastlw_)
       a(ix^s)=wroe(ix^s,mom(idim))-cfast(ix^s)
       jump(ix^s)=half/csound2(ix^s)*(&
            afast(ix^s)*(-cfast(ix^s)*rhodv(ix^s)+dp(ix^s))&
            +aslow(ix^s)*(+cslow(ix^s)*bdv(ix^s)+bdb(ix^s)))
    case(slowrw_)
       a(ix^s)=wroe(ix^s,mom(idim))+cslow(ix^s)
       jump(ix^s)=half/csound2(ix^s)*(&
            aslow(ix^s)*(+cslow(ix^s)*rhodv(ix^s)+dp(ix^s))&
            +afast(ix^s)*(+cfast(ix^s)*bdv(ix^s)-bdb(ix^s)))
    case(slowlw_)
       a(ix^s)=wroe(ix^s,mom(idim))-cslow(ix^s)
       jump(ix^s)=half/csound2(ix^s)*(&
            aslow(ix^s)*(-cslow(ix^s)*rhodv(ix^s)+dp(ix^s))&
            +afast(ix^s)*(-cfast(ix^s)*bdv(ix^s)-bdb(ix^s)))
    case(entrow_)
       a(ix^s)=wroe(ix^s,mom(idim))
       jump(ix^s)=wr(ix^s,rho_)-wl(ix^s,rho_)-dp(ix^s)/csound2(ix^s)
    case(diverw_)
       if(divbwave)then
          a(ix^s)=wroe(ix^s,mom(idim))
          jump(ix^s)=wr(ix^s,mag(idim))-wl(ix^s,mag(idim))
       else
          a(ix^s)=zero
          jump(ix^s)=zero
       endif
    case(alfvrw_)
       a(ix^s)=wroe(ix^s,mom(idim))+dabs(wroe(ix^s,mag(idim)))/wroe(ix^s,rho_)
       jump(ix^s)=+bdv(ix^s)-bdb(ix^s)
    case(alfvlw_)
       a(ix^s)=wroe(ix^s,mom(idim))-dabs(wroe(ix^s,mag(idim)))/wroe(ix^s,rho_)
       jump(ix^s)=-bdv(ix^s)-bdb(ix^s)
    end select
 
    ! Calculate "smalla" or modify "a" based on the "typeentropy" switch
 
    select case(typeentropy(il))
    case('yee')
       ! Based on Yee JCP 68,151 eq 3.23
       smalla(ix^s)=entropycoef(il)
    case('harten','powell', 'ratio')
       ! Based on Harten & Hyman JCP 50, 235 and Zeeuw & Powell JCP 104,56
       ! Initialize left and right eigenvalues by velocities
       al(ix^s)= wl(ix^s,mom(idim))/wl(ix^s,rho_)
       ar(ix^s)= wr(ix^s,mom(idim))/wr(ix^s,rho_)
       ! Calculate the final "aL" and "aR"
       select case(il)
       case(fastrw_)
          ! These quantities will be used for all the fast and slow waves
          ! Calculate soundspeed**2 and cs**2+ca**2.
          call mhd_get_pthermal(wl,x,ixg^ll,ix^l,cs2l)
          cs2l(ix^s)=mhd_gamma*cs2l(ix^s)/wl(ix^s,rho_)
          call mhd_get_pthermal(wr,x,ixg^ll,ix^l,cs2r)
          cs2r(ix^s)=mhd_gamma*cs2r(ix^s)/wl(ix^s,rho_)
          cs2ca2l(ix^s)=cs2l(ix^s)+sum(wl(ix^s,mag(:))**2,dim=ndim+1)/wl(ix^s,rho_)
          cs2ca2r(ix^s)=cs2r(ix^s)+sum(wr(ix^s,mag(:))**2,dim=ndim+1)/wr(ix^s,rho_)
          ! Save the discriminants into cs2L and cs2R
          cs2l(ix^s)=&
               dsqrt(cs2ca2l(ix^s)**2-4d0*cs2l(ix^s)*wl(ix^s,mag(idim))**2/wl(ix^s,rho_))
          cs2r(ix^s)=&
               dsqrt(cs2ca2r(ix^s)**2-4d0*cs2r(ix^s)*wr(ix^s,mag(idim))**2/wr(ix^s,rho_))
 
          ! The left and right eigenvalues for the fast wave going to right
          al(ix^s)=al(ix^s) + dsqrt(half*(cs2ca2l(ix^s) + cs2l(ix^s)))
          ar(ix^s)=ar(ix^s) + dsqrt(half*(cs2ca2r(ix^s) + cs2r(ix^s)))
       case(fastlw_)
          al(ix^s)=al(ix^s) - dsqrt(half*(cs2ca2l(ix^s) + cs2l(ix^s)))
          ar(ix^s)=ar(ix^s) - dsqrt(half*(cs2ca2r(ix^s) + cs2r(ix^s)))
       case(slowrw_)
          al(ix^s)=al(ix^s) + dsqrt(half*(cs2ca2l(ix^s) - cs2l(ix^s)))
          ar(ix^s)=ar(ix^s) + dsqrt(half*(cs2ca2r(ix^s) - cs2r(ix^s)))
       case(slowlw_)
          al(ix^s)=al(ix^s) - dsqrt(half*(cs2ca2l(ix^s) - cs2l(ix^s)))
          ar(ix^s)=ar(ix^s) - dsqrt(half*(cs2ca2r(ix^s) - cs2r(ix^s)))
       case(entrow_,diverw_)
          ! These propagate by the velocity
       case(alfvrw_)
          ! Store the Alfven speeds into cs2ca2L and cs2ca2R
          cs2ca2l(ix^s)=dabs(wl(ix^s,mag(idim)))/dsqrt(wl(ix^s,rho_))
          cs2ca2r(ix^s)=dabs(wr(ix^s,mag(idim)))/dsqrt(wr(ix^s,rho_))
 
          al(ix^s)=al(ix^s) + cs2ca2l(ix^s)
          ar(ix^s)=ar(ix^s) + cs2ca2r(ix^s)
       case(alfvlw_)
          al(ix^s)=al(ix^s) - cs2ca2l(ix^s)
          ar(ix^s)=ar(ix^s) - cs2ca2r(ix^s)
       end select
    end select
 
    call entropyfix(ix^l,il,al,ar,a,smalla)
 
  end subroutine geteigenjump2
 
  ! Multiply q by R(il,iw), where R is the right eigenvalue matrix at wroe
  subroutine mhd_rtimes(q,w,ix^L,iw,il,idim,rq,workroe)
    use mod_global_parameters
 
    integer, intent(in)             :: ix^L, iw, il, idim
    double precision, intent(in)    :: w(ixG^T, nw), q(ixG^T)
    double precision, intent(inout) :: rq(ixG^T)
    double precision, intent(inout) :: workroe(ixG^T, nworkroe)
 
    call rtimes2(q,w,ix^l,iw,il,idim,rq,&
         workroe(ixg^t,1),workroe(ixg^t,2), &
         workroe(ixg^t,3),workroe(ixg^t,4),workroe(ixg^t,5),workroe(ixg^t,6), &
         workroe(ixg^t,7),workroe(ixg^t,14),workroe(ixg^t,15))
 
  end subroutine mhd_rtimes
 
  ! Multiply q by R(il,iw), where R is the right eigenvalue matrix at wroe
  subroutine rtimes2(q,wroe,ix^L,iw,il,idim,rq, &
       cfast,cslow,afast,aslow,csound2,dp,rhodv,bv,v2a2)
    use mod_global_parameters
 
    integer                            :: ix^L,iw,il,idim,idir,jdir
    double precision                   :: wroe(ixG^T,nw)
    double precision, dimension(ixG^T) :: q,rq
    double precision, dimension(ixG^T) :: cfast,cslow,afast,aslow,csound2,dp,rhodv
    double precision, dimension(ixG^T) :: bv,v2a2
 
    idir=idim+1-ndir*(idim/ndir)
    jdir=idir+1-ndir*(idir/ndir)
 
    if(iw == rho_) then
      select case(il)
      case(fastrw_,fastlw_)
        rq(ix^s)=q(ix^s)*afast(ix^s)
      case(slowrw_,slowlw_)
        rq(ix^s)=q(ix^s)*aslow(ix^s)
      case(entrow_)
        rq(ix^s)=q(ix^s)
      case(diverw_,alfvrw_,alfvlw_)
        rq(ix^s)=zero
      end select
    else if(iw == e_) then
      if(il==fastrw_)then
        ! Store 0.5*v**2+(2-gamma)/(gamma-1)*a**2
        v2a2(ix^s)=half*sum(wroe(ix^s,mom(:))**2,dim=ndim+1)+ &
             (two-mhd_gamma)/(mhd_gamma-one)*csound2(ix^s)
        ! Store sgn(bx)*(betay*vy+betaz*vz) in bv
        bv(ix^s)=wroe(ix^s,mag(idir))*wroe(ix^s,mom(idir))
        if(ndir==3)bv(ix^s)=bv(ix^s)+wroe(ix^s,mag(jdir))*wroe(ix^s,mom(jdir))
        bv(ix^s)=bv(ix^s)*sign(one,wroe(ix^s,mag(idim)))
      else if(il==alfvrw_)then
        !Store betaz*vy-betay*vz in bv
        bv(ix^s)=(wroe(ix^s,mag(jdir))*wroe(ix^s,mom(idir))-&
             wroe(ix^s,mag(idir))*wroe(ix^s,mom(jdir)))
      endif
 
      select case(il)
      case(fastrw_)
        rq(ix^s)=q(ix^s)*(-aslow(ix^s)*cslow(ix^s)*bv(ix^s)+afast(ix^s)*&
             (v2a2(ix^s)+cfast(ix^s)*(cfast(ix^s)+wroe(ix^s,mom(idim)))))
      case(fastlw_)
        rq(ix^s)=q(ix^s)*(+aslow(ix^s)*cslow(ix^s)*bv(ix^s)+afast(ix^s)*&
             (v2a2(ix^s)+cfast(ix^s)*(cfast(ix^s)-wroe(ix^s,mom(idim)))))
      case(slowrw_)
        rq(ix^s)=q(ix^s)*(+afast(ix^s)*cfast(ix^s)*bv(ix^s)+aslow(ix^s)*&
             (v2a2(ix^s)+cslow(ix^s)*(cslow(ix^s)+wroe(ix^s,mom(idim)))))
      case(slowlw_)
        rq(ix^s)=q(ix^s)*(-afast(ix^s)*cfast(ix^s)*bv(ix^s)+aslow(ix^s)*&
             (v2a2(ix^s)+cslow(ix^s)*(cslow(ix^s)-wroe(ix^s,mom(idim)))))
      case(entrow_)
        rq(ix^s)= q(ix^s)*half*sum(wroe(ix^s,mom(:))**2,dim=ndim+1)
      case(diverw_)
        if(divbwave)then
          rq(ix^s)= q(ix^s)*wroe(ix^s,mag(idim))
        else
          rq(ix^s)= zero
        endif
      case(alfvrw_)
        rq(ix^s)=+q(ix^s)*bv(ix^s)
      case(alfvlw_)
        rq(ix^s)=-q(ix^s)*bv(ix^s)
      end select
    else if(any(mom(:)==iw)) then
      if(iw==mom(idim))then
        select case(il)
        case(fastrw_)
          rq(ix^s)=q(ix^s)*afast(ix^s)*(wroe(ix^s,iw)+cfast(ix^s))
        case(fastlw_)
          rq(ix^s)=q(ix^s)*afast(ix^s)*(wroe(ix^s,iw)-cfast(ix^s))
        case(slowrw_)
          rq(ix^s)=q(ix^s)*aslow(ix^s)*(wroe(ix^s,iw)+cslow(ix^s))
        case(slowlw_)
          rq(ix^s)=q(ix^s)*aslow(ix^s)*(wroe(ix^s,iw)-cslow(ix^s))
        case(entrow_)
          rq(ix^s)=q(ix^s)*wroe(ix^s,iw)
        case(diverw_,alfvlw_,alfvrw_)
          rq(ix^s)=zero
        end select
      else
        select case(il)
        case(fastrw_)
          rq(ix^s)=q(ix^s)*(afast(ix^s)*wroe(ix^s,iw)-aslow(ix^s)*&
               cslow(ix^s)*wroe(ix^s,mag(1)-mom(1)+iw)*sign(one,wroe(ix^s,mag(idim))))
        case(fastlw_)
          rq(ix^s)=q(ix^s)*(afast(ix^s)*wroe(ix^s,iw)+aslow(ix^s)*&
               cslow(ix^s)*wroe(ix^s,mag(1)-mom(1)+iw)*sign(one,wroe(ix^s,mag(idim))))
        case(slowrw_)
          rq(ix^s)=q(ix^s)*(aslow(ix^s)*wroe(ix^s,iw)+afast(ix^s)*&
               cfast(ix^s)*wroe(ix^s,mag(1)-mom(1)+iw)*sign(one,wroe(ix^s,mag(idim))))
        case(slowlw_)
          rq(ix^s)=q(ix^s)*(aslow(ix^s)*wroe(ix^s,iw)-afast(ix^s)*&
               cfast(ix^s)*wroe(ix^s,mag(1)-mom(1)+iw)*sign(one,wroe(ix^s,mag(idim))))
        case(entrow_)
          rq(ix^s)=q(ix^s)*wroe(ix^s,iw)
        case(diverw_)
          rq(ix^s)=zero
        case(alfvrw_)
          if(iw==mom(idir))then
            rq(ix^s)=+q(ix^s)*wroe(ix^s,mag(jdir))
          else
            rq(ix^s)=-q(ix^s)*wroe(ix^s,mag(idir))
          endif
        case(alfvlw_)
          if(iw==mom(idir))then
            rq(ix^s)=-q(ix^s)*wroe(ix^s,mag(jdir))
          else
            rq(ix^s)=+q(ix^s)*wroe(ix^s,mag(idir))
          endif
        end select
      end if ! iw=m_idir,m_jdir
    else if(any(mag(:)==iw)) then
      if(iw==mag(idim))then
        if(il==diverw_ .and. divbwave)then
          rq(ix^s)=q(ix^s)
        else
          rq(ix^s)=zero
        endif
      else
        select case(il)
        case(fastrw_,fastlw_)
          rq(ix^s)=+q(ix^s)*aslow(ix^s)*dsqrt(csound2(ix^s))*wroe(ix^s,iw)&
               /wroe(ix^s,rho_)
        case(slowrw_,slowlw_)
          rq(ix^s)=-q(ix^s)*afast(ix^s)*dsqrt(csound2(ix^s))*wroe(ix^s,iw)&
               /wroe(ix^s,rho_)
        case(entrow_,diverw_)
          rq(ix^s)=zero
        case(alfvrw_,alfvlw_)
          if(iw==mag(idir))then
            rq(ix^s)=-q(ix^s)*wroe(ix^s,mag(jdir))&
                 /sign(wroe(ix^s,rho_),wroe(ix^s,mag(idim)))
          else
            rq(ix^s)=+q(ix^s)*wroe(ix^s,mag(idir))&
                 /sign(wroe(ix^s,rho_),wroe(ix^s,mag(idim)))
          end if
        end select
      end if ! iw=b_idir,b_jdir
    end if
 
  end subroutine rtimes2
 
end module mod_mhd_roe