C -*- Mode: Fortran; -*-
c-----------------------------------------------------------------------
c     Ravi Samtaney
c     Copyright 2004
c     Princeton Plasma Physics Laboratory
c     All Rights Reserved
c     
c-----------------------------------------------------------------------
c     $Log: inviscidfluxrp.F,v $
c     Revision 1.2  2005/10/11 22:05:30  samtaney
c     Fixed the Roe method.
c
c     Revision 1.1.1.1  2005/09/12 17:28:10  samtaney
c     Original source.
c
c     Revision 1.5  2004/12/04 15:46:16  samtaney
c     Intermediate checkin.
c
c     Revision 1.4  2004/09/24 19:31:59  samtaney
c     Commented the portion of Central Flux computation. Now using upwinding.
c     However, several bugs in the central diff flux were fixed.
c
c     Revision 1.3  2004/09/15 15:45:45  samtaney
c     Uses total pressure now.
c
c     Revision 1.2  2004/06/18 22:30:14  samtaney
c     Intermediate checkin.
c
c     Revision 1.1.1.1  2004/05/28 19:33:27  samtaney
c     Original source - copied from 3D Unsplit code.
c
c     Revision 1.2  2004/04/12 18:08:34  samtaney
c     Normal component in RP solution is average of left/right states.
c
c-----------------------------------------------------------------------
c     Computes fluxes using the ROE method. 
c-----------------------------------------------------------------------
      subroutine InviscidFluxXi(finv,vx)
      use mesh
      use mesh_common
      use properties
      use MagAxisToroidalField

      implicit none
c     
c     This subroutine determines the vector of inviscid fluxes.
      double precision:: finv(1:nxlocal+1,1:nylocal,1:nzlocal,nvar)
      double precision,intent(in):: 
     &     vx(IXLO:IXHI,IYLO:IYHI,IZLO:IZHI,nvar)
c      double precision,intent(in):: 
c     &     beq(IXLO:IXHI,IYLO:IYHI,IZLO:IZHI,3)
      integer:: direction,ilo,ihi
      integer:: i,j,k,l
c     
      double precision:: vd(IXLO:IXHI,nvar)
      double precision:: vl(1:nxlocal+1,nvar)
      double precision:: vr(1:nxlocal+1,nvar)
      double precision:: wl(1:nxlocal+1,nvar)
      double precision:: wr(1:nxlocal+1,nvar)
      double precision:: btt(1:nxlocal+1)
      double precision:: fn(1:nxlocal+1,nvar)
      double precision:: fxi(1:nxlocal+1,nvar)
      double precision:: fc(1:nxlocal+1,nvar)
      double precision:: vrp(1:nxlocal+1,nvar)
      double precision:: Rclocal(0:nxlocal+1)
      double precision:: costheta(1:nxlocal+1)
      double precision:: sintheta(1:nxlocal+1)
      double precision:: costhetaC(0:nxlocal+1)
      double precision:: sinthetaC(0:nxlocal+1)
      double precision:: dlen(1:nxlocal+1)
c     
      ilo=1+hasxlo
      ihi=nxlocal+1-hasxhi
c      write(6,*) 'INVISCIDFLUX XI',ilo,ihi, hasxlo
      do k=1,nzlocal,1
         do j=1,nylocal,1
            vd(ixlo:ixhi,:)=vx(ixlo:ixhi,j,k,:)
            btt(ilo:ihi)=Btaxis/RXi(ilo:ihi,j)
            costheta(ilo:ihi)=dZdEta(ilo:ihi,j)
     &           /dlXi(ilo:ihi,j)
            sintheta(ilo:ihi)=-dRdEta(ilo:ihi,j)
     &           /dlXi(ilo:ihi,j)
            call ConstructLRStatesXi(vd,vl,vr)
            call RemapVariablesToLocalXi(vl,wl,costheta, sintheta,
     &           ilo,ihi)
            call RemapVariablesToLocalXi(vr,wr,costheta, sintheta,
     &           ilo,ihi)
            call InviscidFluxRP(fn,vrp,wl,wr,btt,1,nxlocal+1,
     &           ilo,ihi,nvar,1)
c
            call RemapFluxXi(fn,fxi,costheta,sintheta,ilo,ihi)
            do l=1,nvar,1
               finv(ilo:ihi,j,k,l)=RXi(ilo:ihi,j)*fxi(ilo:ihi,l)
     &              *dlXi(ilo:ihi,j)
            enddo
            l=7
            finv(ilo:ihi,j,k,l)=fxi(ilo:ihi,l)
     &              *dlXi(ilo:ihi,j)
c
         enddo
      enddo
c     
      return
      end
c
c-----------------------------------------------------------------------
         subroutine RemapVariablesToLocalXi(v,w,costheta,sintheta,
     &     ilo,ihi)
         use mesh_parms
         implicit none
         integer:: ilo,ihi
         double precision:: v(1:nxlocal+1,nvar)
         double precision:: w(1:nxlocal+1,nvar)
         double precision:: costheta(1:nxlocal+1)
         double precision:: sintheta(1:nxlocal+1)
c
         w(ilo:ihi,1)=v(ilo:ihi,1)
         w(ilo:ihi,2)=(v(ilo:ihi,2)*costheta(ilo:ihi)+
     &        v(ilo:ihi,3)*sintheta(ilo:ihi))
         w(ilo:ihi,3)=(-v(ilo:ihi,2)*sintheta(ilo:ihi)+
     &        v(ilo:ihi,3)*costheta(ilo:ihi))
         w(ilo:ihi,4)=v(ilo:ihi,4)
c         w(ilo:ihi,4)=-v(ilo:ihi,4)
         w(ilo:ihi,5)=(v(ilo:ihi,5)*costheta(ilo:ihi)+
     &        v(ilo:ihi,6)*sintheta(ilo:ihi))
c     KLUDGE
c         w(ilo:ihi,5)=0.D0
         w(ilo:ihi,6)=(-v(ilo:ihi,5)*sintheta(ilo:ihi)+
     &        v(ilo:ihi,6)*costheta(ilo:ihi))
         w(ilo:ihi,7:8)=v(ilo:ihi,7:8)
         return
         end
c-----------------------------------------------------------------------
         subroutine RemapFluxXi(fn,fxi,costheta,sintheta,
     &     ilo,ihi)
         use mesh_parms
         implicit none
         integer:: ilo,ihi
         double precision:: fn(1:nxlocal+1,nvar)
         double precision:: fxi(1:nxlocal+1,nvar)
         double precision:: costheta(1:nxlocal+1)
         double precision:: sintheta(1:nxlocal+1)
c
         fxi(ilo:ihi,1)=fn(ilo:ihi,1)
         fxi(ilo:ihi,2)=fn(ilo:ihi,2)*costheta(ilo:ihi)-
     &        fn(ilo:ihi,3)*sintheta(ilo:ihi)
         fxi(ilo:ihi,3)=fn(ilo:ihi,2)*sintheta(ilo:ihi)+
     &        fn(ilo:ihi,3)*costheta(ilo:ihi)
         fxi(ilo:ihi,4)=fn(ilo:ihi,4)
         fxi(ilo:ihi,5)=fn(ilo:ihi,5)*costheta(ilo:ihi)-
     &        fn(ilo:ihi,6)*sintheta(ilo:ihi)
         fxi(ilo:ihi,6)=fn(ilo:ihi,5)*sintheta(ilo:ihi)+
     &        fn(ilo:ihi,6)*costheta(ilo:ihi)
         fxi(ilo:ihi,7:8)=fn(ilo:ihi,7:8)
c
         return
         end
        
c-----------------------------------------------------------------------
      subroutine InviscidFluxEta(finv,vx)
      use mesh
      use mesh_common
      use properties
      use MagAxisToroidalField

      implicit none
c     
c     This subroutine determines the vector of inviscid fluxes.
      double precision:: finv(1:nxlocal,1:nylocal+1,1:nzlocal,nvar)
      double precision,intent(in):: 
     &     vx(IXLO:IXHI,IYLO:IYHI,IZLO:IZHI,nvar)
      integer:: direction,ilo,ihi
      integer:: i,j,k,l
c     
      double precision:: vd(IYLO:IYHI,nvar)
      double precision:: vl(1:nylocal+1,nvar)
      double precision:: vr(1:nylocal+1,nvar)
      double precision:: wl(1:nylocal+1,nvar)
      double precision:: wr(1:nylocal+1,nvar)
      double precision:: fn(1:nylocal+1,nvar)
      double precision:: feta(1:nylocal+1,nvar)
      double precision:: fc(1:nylocal+1,nvar)
      double precision:: vrp(1:nylocal+1,nvar)
      double precision:: Rclocal(0:nylocal+1)
      double precision:: costheta(1:nylocal+1)
      double precision:: sintheta(1:nylocal+1)
      double precision:: costhetaC(0:nylocal+1)
      double precision:: sinthetaC(0:nylocal+1)
      double precision:: btt(1:nylocal+1)
c     
      ilo=1+hasylo
      ihi=nylocal+1-hasyhi
c      write(6,*) 'INVISCIDFLUX ETA',ilo,ihi
      do k=1,nzlocal,1
         do i=1,nxlocal,1
            vd(iylo:iyhi,:)=vx(i,iylo:iyhi,k,:)
            btt(ilo:ihi)=Btaxis/REta(i,ilo:ihi)
            costheta(ilo:ihi)=-dZdXi(i,ilo:ihi)
     &           /dlEta(i,ilo:ihi)
            sintheta(ilo:ihi)=dRdXi(i,ilo:ihi)
     &           /dlEta(i,ilo:ihi)
            call ConstructLRStatesEta(vd,vl,vr)
            call RemapVariablesToLocalEta(vl,wl,costheta, sintheta,
     &           ilo,ihi)
            call RemapVariablesToLocalEta(vr,wr,costheta, sintheta,
     &           ilo,ihi)
            call InviscidFluxRP(fn,vrp,wl,wr,btt,1,nylocal+1,
     &           ilo,ihi,nvar,2)
c
            call RemapFluxEta(fn,feta,costheta,sintheta,ilo,ihi)
            do l=1,nvar,1
               finv(i,ilo:ihi,k,l)=REta(i,ilo:ihi)*feta(ilo:ihi,l)
     &              *dlEta(i,ilo:ihi)
            enddo
            l=7
            finv(i,ilo:ihi,k,l)=feta(ilo:ihi,l)
     &           *dlEta(i,ilo:ihi)
         enddo
      enddo
c     
      return
      end
c
c-----------------------------------------------------------------------
         subroutine RemapVariablesToLocalEta(v,w,costheta,sintheta
     &     ,ilo,ihi)
         use mesh_parms
         implicit none
         integer:: ilo,ihi
         double precision:: v(1:nylocal+1,nvar)
         double precision:: w(1:nylocal+1,nvar)
         double precision:: costheta(1:nylocal+1)
         double precision:: sintheta(1:nylocal+1)
c
         w(ilo:ihi,1)=v(ilo:ihi,1)
         w(ilo:ihi,2)=v(ilo:ihi,2)*costheta(ilo:ihi)+
     &        v(ilo:ihi,3)*sintheta(ilo:ihi)
         w(ilo:ihi,3)=-v(ilo:ihi,2)*sintheta(ilo:ihi)+
     &        v(ilo:ihi,3)*costheta(ilo:ihi)
         w(ilo:ihi,4)=v(ilo:ihi,4)
c         w(ilo:ihi,4)=-v(ilo:ihi,4)
         w(ilo:ihi,5)=v(ilo:ihi,5)*costheta(ilo:ihi)+
     &        v(ilo:ihi,6)*sintheta(ilo:ihi)
         w(ilo:ihi,6)=-v(ilo:ihi,5)*sintheta(ilo:ihi)+
     &        v(ilo:ihi,6)*costheta(ilo:ihi)
cKLUDGE
c         w(ilo:ihi,6)=0.D0
         w(ilo:ihi,7)=v(ilo:ihi,7)
         w(ilo:ihi,8)=v(ilo:ihi,8)
         return
         end
c-----------------------------------------------------------------------
         subroutine RemapFluxEta(fn,feta,costheta,sintheta,ilo,ihi)
         use mesh
         integer:: ilo,ihi
         double precision:: fn(1:nylocal+1,nvar)
         double precision:: feta(1:nylocal+1,nvar)
         double precision:: costheta(1:nylocal+1)
         double precision:: sintheta(1:nylocal+1)
c
         feta(ilo:ihi,1)=fn(ilo:ihi,1)
         feta(ilo:ihi,2)=fn(ilo:ihi,2)*costheta(ilo:ihi)-
     &        fn(ilo:ihi,3)*sintheta(ilo:ihi)
         feta(ilo:ihi,3)=fn(ilo:ihi,2)*sintheta(ilo:ihi)+
     &        fn(ilo:ihi,3)*costheta(ilo:ihi)
         feta(ilo:ihi,4)=fn(ilo:ihi,4)
         feta(ilo:ihi,5)=fn(ilo:ihi,5)*costheta(ilo:ihi)-
     &        fn(ilo:ihi,6)*sintheta(ilo:ihi)
         feta(ilo:ihi,6)=fn(ilo:ihi,5)*sintheta(ilo:ihi)+
     &        fn(ilo:ihi,6)*costheta(ilo:ihi)
         feta(ilo:ihi,7:8)=fn(ilo:ihi,7:8)
c         feta(ilo:ihi,8)=fn(ilo:ihi,8)
         return
         end
c-----------------------------------------------------------------------
      subroutine InviscidFluxPhi(finv,vx)
      use mesh
      use mesh_common
      use properties
      use MagAxisToroidalField
c
      implicit none

c     
c     This subroutine determines the vector of inviscid fluxes.
      double precision:: finv(1:nxlocal,1:nylocal,1:nzlocal+1,nvar)
      double precision,intent(in):: 
     &     vx(IXLO:IXHI,IYLO:IYHI,IZLO:IZHI,nvar)
      integer:: direction,ilo,ihi
      integer:: i,j,k,l
c     
      double precision:: vd(IZLO:IZHI,nvar)
      double precision:: vl(1:nzlocal+1,nvar)
      double precision:: vr(1:nzlocal+1,nvar)
      double precision:: wl(1:nzlocal+1,nvar)
      double precision:: wr(1:nzlocal+1,nvar)
      double precision:: fn(1:nzlocal+1,nvar)
      double precision:: fphi(1:nzlocal+1,nvar)
      double precision:: fc(1:nzlocal+1,nvar)
      double precision:: vrp(1:nzlocal+1,nvar)
      double precision:: btt(1:nzlocal+1)
c     
      ilo=1
      ihi=nzlocal+1
c      write(6,*) 'INVISCIDFLUX PHI',ilo,ihi
      do j=1,nylocal,1
         do i=1,nxlocal,1
            vd(izlo:izhi,:)=vx(i,j,izlo:izhi,:)
            btt(ilo:ihi)=Btaxis/Rc(i,j)
            call ConstructLRStatesPhi(vd,vl,vr)
            call RemapVariablesToLocalPhi(vl,wl,ilo,ihi)
            call RemapVariablesToLocalPhi(vr,wr,ilo,ihi)
            call InviscidFluxRP(fn,vrp,wl,wr,btt,1,nzlocal+1,
     &           ilo,ihi,nvar,3)
            call RemapFluxPhi(fn,fphi,ilo,ihi)
            finv(i,j,ilo:ihi,:)=fphi(ilo:ihi,:)
         enddo
      enddo
c     
      return
      end
c
c-----------------------------------------------------------------------
         subroutine RemapVariablesToLocalPhi(v,w,ilo,ihi)
         use mesh
         implicit none
         integer:: ind
         integer:: ilo,ihi
         double precision:: v(1:nzlocal+1,nvar)
         double precision:: w(1:nzlocal+1,nvar)
c
         w(ilo:ihi,1)=v(ilo:ihi,1)
         w(ilo:ihi,2)=v(ilo:ihi,4)
c         w(ilo:ihi,2)=-v(ilo:ihi,4)
         w(ilo:ihi,3)=v(ilo:ihi,2)
         w(ilo:ihi,4)=v(ilo:ihi,3)
         w(ilo:ihi,5)=v(ilo:ihi,7)
         w(ilo:ihi,6)=v(ilo:ihi,5)
         w(ilo:ihi,7)=v(ilo:ihi,6)
         w(ilo:ihi,8)=v(ilo:ihi,8)
         return
         end
c-----------------------------------------------------------------------
         subroutine RemapFluxPhi(fn,fphi,ilo,ihi)
         use mesh
         implicit none
         integer:: ilo,ihi
         double precision:: fn(1:nzlocal+1,nvar)
         double precision:: fphi(1:nzlocal+1,nvar)
c
         fphi(ilo:ihi,1)=fn(ilo:ihi,1)
         fphi(ilo:ihi,4)=fn(ilo:ihi,2)
         fphi(ilo:ihi,2)=fn(ilo:ihi,3)
         fphi(ilo:ihi,3)=fn(ilo:ihi,4)
         fphi(ilo:ihi,7)=fn(ilo:ihi,5)
         fphi(ilo:ihi,5)=fn(ilo:ihi,6)
         fphi(ilo:ihi,6)=fn(ilo:ihi,7)
         fphi(ilo:ihi,8)=fn(ilo:ihi,8)
         return
         end
c-----------------------------------------------------------------------
      subroutine InviscidFluxRP(finv,vm,vl,vr,btt,inlo,
     &     inhi,ilo,ihi,nvar,direction)
c      use mesh
c      use mesh_common
      use properties
      use MagAxisToroidalField
      use alphas
      use evectors
      double precision:: finv(inlo:inhi,nvar)
      double precision,intent(in):: vl(inlo:inhi,nvar)
      double precision,intent(in):: vr(inlo:inhi,nvar)
      double precision,intent(in):: btt(inlo:inhi)
      double precision:: vm(inlo:inhi,nvar)
      double precision:: rho,u,v,w,bx,by,bz
      double precision:: eint, lam, press
        double precision:: alamda(ILO:IHI,nvar)
	double precision:: evl(ILO:IHI,nvar,nvar)
	double precision:: evr(ILO:IHI,nvar,nvar)
	double precision:: durl(ILO:IHI,nvar)
	double precision:: alpha(ILO:IHI,nvar)
	double precision:: bxsq,bysq,bzsq
        integer:: i,j,k,l,m,ilo,ihi,direction,ii
        double precision:: Btlocal
c
c	 This subroutine determines the vector of inviscid fluxes.
c
c	 v: Vector of variables - rho, u,v,w,Bx,By,Bz,p
c	    rho=density, u=x-velocity, v=y-velocity
c           w=z-velocity
c           Bx,By,Bz= Magnetic Field
c	    p=pressure
c
c	 vl(i) = Left state after linear representation in cell, i .
c		 at the i-th interface
c	 vr(i) = Right state after linear representation in cell, i+1 .
c		 at the i-th interface
c
c
      double precision:: utilde(ilo:ihi,nvar)
c	 utilde = vector of Roe averaged variables.
c	 utilde(1,2,3,4,5,6,7,8) = {rhotilde, utilde, vtilde, wtilde, 
c                           Bxtilde,Bytilde,Bztilde,ptilde}
c
c	 Set the Roe averaged variables depending  on vl,vr
c      write(6,*) 'In Inviscidflux Roe', indx,indy,indz
	 call SetRoeVariables(utilde,vl,vr,inlo,inhi,ilo,ihi)
c
c         if(direction.eq.2) then
            Btlocal=0.D0
c         else
c            Btlocal=Btaxis
c         endif
c
        alamda=0.D0
        durl=0.D0
c
	 call SetEigenSystem(utilde,evl,evr,alamda,ilo,ihi,
     &       btt,inlo,inhi,direction,nvar)
c	 call SetDurl(durl,vl,vr,ilo,ihi)
         durl(ilo:ihi,:)=vr(ilo:ihi,:)-vl(ilo:ihi,:)
	 call SetAlphas(utilde,alpha,durl,evl,evr,ilo,ihi)
         call RiemannProblem(evl,evr,alamda,alpha,durl,vl,vr,vm,
     &        inlo,inhi,ilo,ihi)
c
c         vm=0.5D0*(vl+vr)
         do i=ilo,ihi,1
            rho=vm(i,1)
            u=vm(i,2); v=vm(i,3); w=vm(i,4)
            bx=vm(i,5); by=vm(i,6); bz=vm(i,7)
            bxsq=bx*bx
            bysq=by*by
            bzsq=bz*bz
            press=vm(i,8)-0.5D0*(bxsq+bysq+bzsq)
            finv(i,1)=rho*u
            finv(i,2)=rho*u*u +press+0.5D0*(bysq+bzsq-bxsq)
c            finv(i,2)=rho*u*u -bxsq
c     &           Btaxis*bz
            finv(i,3)=rho*u*v-bx*by
            finv(i,4)=rho*u*w-bx*bz
c            finv(i,4)=rho*u*w-bx*(bz+Btlocal)
c            finv(i,2)=press+0.5D0*(bysq+bzsq-bxsq)
c            finv(i,3)=-bx*by
c            finv(i,4)=-bx*bz
c            eint=vm(i,8)/(gamma-1.D0)
            eint=press/(gamma-1.D0)
            finv(i,5)=0.D0
            finv(i,6)=u*by-v*bx
c            finv(i,7)=u*(bz+Btaxis)-w*bx
            finv(i,7)=u*(bz)-w*bx
            finv(i,8)=(0.5D0*rho*(u*u+v*v+w*w)+eint+
     &           press+(bxsq+bysq+bzsq))*vm(i,2)-
     &           bx*(u*bx+v*by+w*bz)
c     &           press+(bxsq+bysq+bzsq)+Btaxis*bz)*vm(i,2)-
c$$$            finv(i,8)=(0.5D0*rho*(u*u+v*v+w*w)+eint+
c$$$     &           press+(bxsq+bysq+bzsq))*vm(i,2)
c
            if(direction.eq.1.or.direction.eq.2) then
               finv(i,2)=finv(i,2)+Btt(i)*bz
               finv(i,7)=finv(i,7)+Btt(i)*u
               finv(i,8)=finv(i,8)+Btt(i)*bz*u
            else if(direction.eq.3) then
               finv(i,3)=finv(i,3)-Btt(i)*by
               finv(i,4)=finv(i,4)-Btt(i)*bz
               finv(i,6)=finv(i,6)-Btt(i)*v
               finv(i,7)=finv(i,7)-Btt(i)*w
               finv(i,8)=finv(i,8)-Btt(i)*(by*v+bz*w)
            endif
         enddo
c
c
	return
        end

c-----------------------------------------------------------------------
      subroutine RiemannProblem(evl,evr,alamda,alpha,durl,vl,vr,vm,
     &    inlo,inhi,ilo,ihi)

      use mesh
      use mesh_common
      use properties
      integer:: inlo,inhi,ilo,ihi
      integer:: i,j,k,l,m,direction
      double precision:: vl(inlo:inhi,nvar)
      double precision:: vr(inlo:inhi,nvar)
      double precision:: vml(inlo:inhi,nvar)
      double precision:: vmr(inlo:inhi,nvar)
      double precision:: vm(inlo:inhi,nvar)
      double precision:: durl(ilo:ihi,nvar)
      double precision:: rho,u,v,w,bx,by,bz
      double precision:: eint, lam
        double precision:: alamda(ILO:IHI,nvar)
	double precision:: evl(ILO:IHI,nvar,nvar)
	double precision:: evr(ILO:IHI,nvar,nvar)
	double precision:: alpha(ILO:IHI,nvar)
	double precision:: tmp(nvar,nvar)
c
      vmr(ilo:ihi,:)=vr(ilo:ihi,:)
      vml(ilo:ihi,:)=vl(ilo:ihi,:)
      do i=ilo,ihi,1
         do k=1,nvar,1
            if(alamda(i,k).gt.0.D0) then
               vmr(i,:)=vmr(i,:)-alpha(i,k)*evr(i,:,k)
            else if(alamda(i,k).lt.0.D0) then
               vml(i,:)=vml(i,:)+alpha(i,k)*evr(i,:,k)
            else
               vmr(i,:)=vmr(i,:)-alpha(i,k)*evr(i,:,k)
               vml(i,:)=vml(i,:)+alpha(i,k)*evr(i,:,k)
            endif
         enddo
      enddo
c
       vm(ilo:ihi,:)=0.5D0*(vml(ilo:ihi,:)+vmr(ilo:ihi,:))
       vm(ilo:ihi,5)=0.5D0*(vl(ilo:ihi,5)+vr(ilo:ihi,5))
C      vm(ilo:ihi,:)=0.5D0*(vl(ilo:ihi,:)+vr(ilo:ihi,:))
      return
      end
c-----------------------------------------------------------------------
c      subroutine InviscidFluxRoe(finv,vl,vr,ilo,ihi,direction)
c$$$      use mesh
c$$$      use mesh_common
c$$$      use properties
c$$$      use slamdas
c$$$      double precision:: finv(inlo:inhi,nvar)
c$$$      double precision:: vl(inlo:inhi,nvar)
c$$$      double precision:: vr(inlo:inhi,nvar)
c$$$      double precision:: sls(inlo:inhi,nvar)
c$$$      double precision:: finvl(nvar), finvr(nvar)
c$$$      double precision:: rho,u,v,w,bx,by,bz
c$$$      double precision:: eint, lam
c$$$      integer:: direction,ilo,ihi
c$$$      integer:: i,j,k,l
c$$$c	 This subroutine determines the vector of inviscid fluxes.
c$$$c
c$$$c	 v: Vector of variables - rho, u,v,w,Bx,By,Bz,p
c$$$c	    rho=density, u=x-velocity, v=y-velocity
c$$$c           w=z-velocity
c$$$c           Bx,By,Bz= Magnetic Field
c$$$c	    p=pressure
c$$$c
c$$$c	 vl(i) = Left state after linear representation in cell, i .
c$$$c		 at the i-th interface
c$$$c	 vr(i) = Right state after linear representation in cell, i+1 .
c$$$c		 at the i-th interface
c$$$c
c$$$c
c$$$      double precision:: utilde(inlo:inhi,nvar)
c$$$c	 utilde = vector of Roe averaged variables.
c$$$c	 utilde(1,2,3,4,5,6,7,8) = {rhotilde, utilde, vtilde, wtilde, 
c$$$c                           Bxtilde,Bytilde,Bztilde,ptilde}
c$$$c	 sls = vector of S |Lamda| S^-1
c$$$c	 S: matrix of left eigenvectors
c$$$c	 |Lamda|: absolute value of eigenvalues
c$$$c	 S^-1: Inverse of S.
c$$$c	 
c$$$c	 Fluxes finvl=f(u_l)
c$$$c		finvr=f(u_r)
c$$$c
c$$$c
c$$$c	 Set the Roe averaged variables depending  on vl,vr
c$$$c      write(6,*) 'In Inviscidflux Roe', indx,indy,indz
c$$$	 call SetRoeVariables(utilde,vl,vr,ilo,ihi)
c$$$c
c$$$	 call SetSLS(utilde,sls,vl,vr,ilo,ihi,direction)
c$$$c
c$$$         do i=ilo,ihi,1
c$$$            rho=vl(i,1)
c$$$            u=vl(i,2); v=vl(i,3); w=vl(i,4)
c$$$            bx=vl(i,5); by=vl(i,6); bz=vl(i,7)
c$$$            finvl(1)=rho*u
c$$$            finvl(2)=rho*u*u +vl(i,8)+0.5D0*(by*by+bz*bz-bx*bx)
c$$$            finvl(3)=rho*u*v-bx*by
c$$$            finvl(4)=rho*u*w-bx*bz
c$$$            eint=vl(i,8)/(gamma-1.D0)
c$$$            finvl(5)=0.D0
c$$$            finvl(6)=u*by-v*bx
c$$$            finvl(7)=u*bz-w*bx
c$$$            finvl(8)=(0.5D0*rho*(u*u+v*v+w*w)+eint+
c$$$     &           vl(i,8)+(bx*bx+by*by+bz*bz))*vl(i,2)-
c$$$     &           bx*(u*bx+v*by+w*bz)
c$$$c
c$$$            rho=vr(i,1)
c$$$            u=vr(i,2); v=vr(i,3); w=vr(i,4)
c$$$            bx=vr(i,5); by=vr(i,6); bz=vr(i,7)
c$$$            finvr(1)=rho*u
c$$$            finvr(2)=rho*u*u +vr(i,8)+0.5D0*(by*by+bz*bz-bx*bx)
c$$$            finvr(3)=rho*u*v-bx*by
c$$$            finvr(4)=rho*u*w-bx*bz
c$$$            eint=vr(i,8)/(gamma-1.D0)
c$$$            finvr(5)=0.D0
c$$$            finvr(6)=u*by-v*bx
c$$$            finvr(7)=u*bz-w*bx
c$$$            finvr(8)=(0.5D0*rho*(u*u+v*v+w*w)+eint+
c$$$     &           vr(i,8)+(bx*bx+by*by+bz*bz))*vr(i,2)-
c$$$     &           bx*(u*bx+v*by+w*bz)
c$$$c
c$$$            finv(i,1:nvar)=0.5D0*(finvl(1:nvar)+finvr(1:nvar))
c$$$     &           -0.5D0*sls(i,1:nvar)
c$$$         enddo
c$$$c
c$$$c
c$$$	return
c$$$        end
c$$$