C -*- Mode: Fortran; -*-
c-----------------------------------------------------------------------
c     Ravi Samtaney
c     Copyright 2004
c     Princeton Plasma Physics Laboratory
c     All Rights Reserved
c     
c-----------------------------------------------------------------------
c     $Log$
c-----------------------------------------------------------------------
c     Computes fluxes using the HLL 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,nvar)
      double precision,intent(in):: 
     &     vx(IXLO:IXHI,IYLO:IYHI,nvar)
c      double precision,intent(in):: 
c     &     beq(IXLO:IXHI,IYLO:IYHI,IZLO:IZHI,3)
      integer:: direction,ilo,ihi
      integer:: i,j,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 j=1,nylocal,1
            vd(ixlo:ixhi,:)=vx(ixlo:ixhi,j,:)
            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,l)=RXi(ilo:ihi,j)*fxi(ilo:ihi,l)
     &              *dlXi(ilo:ihi,j)
c$$$               do i=ilo,ihi,1
c$$$                  write(6,*) 'FXI=',i,j,l,finv(i,j,l)
c$$$               enddo
            enddo
c
         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)*costheta(ilo:ihi)+
     &        v(ilo:ihi,2)*sintheta(ilo:ihi))
         w(ilo:ihi,2)=(-v(ilo:ihi,1)*sintheta(ilo:ihi)+
     &        v(ilo:ihi,2)*costheta(ilo:ihi))
         w(:,1)=0.D0
         w(ilo:ihi,3:4)=v(ilo:ihi,3:4)
         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)*costheta(ilo:ihi)-
     &        fn(ilo:ihi,2)*sintheta(ilo:ihi)
         fxi(ilo:ihi,2)=fn(ilo:ihi,1)*sintheta(ilo:ihi)+
     &        fn(ilo:ihi,2)*costheta(ilo:ihi)
         fxi(ilo:ihi,3:4)=fn(ilo:ihi,3:4)
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,nvar)
      double precision,intent(in):: 
     &     vx(IXLO:IXHI,IYLO:IYHI,nvar)
      integer:: direction,ilo,ihi
      integer:: i,j,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 i=1,nxlocal,1
            vd(iylo:iyhi,:)=vx(i,iylo:iyhi,:)
            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,l)=REta(i,ilo:ihi)*feta(ilo:ihi,l)
     &              *dlEta(i,ilo:ihi)
c$$$               do j=ilo,ihi,1
c$$$                  write(6,*) 'FETA=',i,j,l,finv(i,j,l)
c$$$               enddo
            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)*costheta(ilo:ihi)+
     &        v(ilo:ihi,2)*sintheta(ilo:ihi)
         w(ilo:ihi,2)=-v(ilo:ihi,1)*sintheta(ilo:ihi)+
     &        v(ilo:ihi,2)*costheta(ilo:ihi)
         w(ilo:ihi,2)=0.D0
         w(ilo:ihi,3:4)=v(ilo:ihi,3:4)
         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)*costheta(ilo:ihi)-
     &        fn(ilo:ihi,2)*sintheta(ilo:ihi)
         feta(ilo:ihi,2)=fn(ilo:ihi,1)*sintheta(ilo:ihi)+
     &        fn(ilo:ihi,2)*costheta(ilo:ihi)
         feta(ilo:ihi,3:4)=fn(ilo:ihi,3:4)
         return
         end
c-----------------------------------------------------------------------
      subroutine InviscidFluxRP(finv,vm,vl,vr,btt,inlo,
     &     inhi,ilo,ihi,nvar,direction)
      use properties
      use MagAxisToroidalField
c      use alphas
c      use evalues
      implicit none
      integer:: inlo,inhi, ilo,ihi,nvar
      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:: alamdaL(ILO:IHI,2)
        double precision:: alamdaR(ILO:IHI,2)
	double precision:: durl(ILO:IHI,nvar)
	double precision:: alpha(ILO:IHI,nvar)
        double precision:: finvl(nvar)
        double precision:: finvr(nvar)
	double precision:: bxsq,bysq,bzsq
        integer:: i,j,l,m,direction,ii
        double precision:: Btlocal
        double precision:: lamdaMax
        double precision:: lamdaMin
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=total pressure 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	 call SetRoeVariables(utilde,vl,vr,inlo,inhi,ilo,ihi)
c
c        alamdaL=0.D0
c        alamdaR=0.D0
        durl=0.D0
c
c$$$        utilde(ilo:ihi,:)=vl(ilo:ihi,:)
c$$$	 call SetEigenValues(utilde,alamdaL,ilo,ihi,
c$$$     &       btt,inlo,inhi,direction,nvar)
c$$$        utilde(ilo:ihi,:)=vr(ilo:ihi,:)
c$$$	 call SetEigenValues(utilde,alamdaR,ilo,ihi,
c$$$     &       btt,inlo,inhi,direction,nvar)
c$$$	 call SetDurl(durl,vl,vr,inlo,inhi,ilo,ihi)
c
         do i=ilo,ihi,1
c            lamdaMin=dmin1(alamdaL(i,2),alamdaR(i,2))
c            lamdaMax=dmax1(alamdaL(i,1),alamdaR(i,1))
c
            bx=vl(i,1); by=vl(i,2); bz=vl(i,3)
            bxsq=bx*bx
            bysq=by*by
            bzsq=bz*bz
            press=vl(i,4)-0.5D0*(bysq+bzsq+bxsq)
            finvl(1)=press+0.5D0*(bysq+bzsq-bxsq)
            finvl(2)=-bx*by
            finvl(3)=-bx*bz
            finvl(4)=0.D0
c
            if(direction.eq.1.or.direction.eq.2) then
               finvl(1)=finvl(1)+Btt(i)*bz
            endif
c
            bx=vr(i,1); by=vr(i,2); bz=vr(i,3)
            bxsq=bx*bx
            bysq=by*by
            bzsq=bz*bz
            press=vr(i,4)-0.5D0*(bysq+bzsq+bxsq)
            finvr(1)=press+0.5D0*(bysq+bzsq-bxsq)
            finvr(2)=-bx*by
            finvr(3)=-bx*bz
            finvr(4)=0.D0
c     
            if(direction.eq.1.or.direction.eq.2) then
               finvr(1)=finvr(1)+Btt(i)*bz
            endif
c$$$            if(lamdaMin.gt.0.D0) then
c$$$               finv(i,1:8)=finvl(1:8)
c$$$            else if(lamdaMax.lt.0.D0) then
c$$$               finv(i,1:8)=finvr(1:8)
c$$$            else if(lamdaMin.ne.lamdaMax) then
c$$$               finv(i,1:8)=(lamdaMax*finvl(1:8)-lamdaMin*finvr(1:8)
c$$$     &                     +lamdaMin*lamdaMax*durl(i,1:8))/
c$$$     &              (lamdaMax-lamdaMin)
c$$$            else
               finv(i,1:4)=0.5D0*(finvl(1:4)+finvr(1:4))
c$$$            endif
         enddo
c
c
	return
        end