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 
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(0: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)
      double precision:: finvl(nvar), finvr(nvar)
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,:)
c            btt(ilo:ihi)=Btaxis/RXi(ilo:ihi,j)
            btt(ilo-1:ihi)=Btaxis/RC(ilo-1:ihi,j)
            call ConstructLRStatesXi(vd,vl,vr)
c
         do i=ilo,ihi,1
c
c$$$            finvl(1)=Rxi(i,j)*dZdEta(i,j)*(vl(i,4)+Btt(i)*vl(i,3))
c$$$            finvl(2)=-Rxi(i,j)*dRdEta(i,j)*(vl(i,4)+Btt(i)*vl(i,3))
c$$$            finvl(3)=0.D0
c$$$            finvl(4)=0.D0
c$$$c
c$$$            finvr(1)=Rxi(i,j)*dZdEta(i,j)*(vr(i,4)+Btt(i)*vr(i,3))
c$$$            finvr(2)=-Rxi(i,j)*dRdEta(i,j)*(vr(i,4)+Btt(i)*vr(i,3))
c$$$            finvr(3)=0.D0
c$$$            finvr(4)=0.D0
c
            finvl(1)=RC(i-1,j)*dZdEtaC(i-1,j)*(vl(i,4)+Btt(i-1)*vl(i,3))
            finvl(2)=-RC(i-1,j)*dRdEtaC(i-1,j)*
     &           (vl(i,4)+Btt(i-1)*vl(i,3))
            finvl(3)=0.D0
            finvl(4)=0.D0
c
            finvr(1)=RC(i,j)*dZdEtaC(i,j)*(vr(i,4)+Btt(i)*vr(i,3))
            finvr(2)=-RC(i,j)*dRdEtaC(i,j)*(vr(i,4)+Btt(i)*vr(i,3))
            finvr(3)=0.D0
            finvr(4)=0.D0
c     
            finv(i,j,1:4)=0.5D0*(finvl(1:4)+finvr(1:4))
         enddo
      enddo
c
c     
      return
      end
c
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(0:nylocal+1)
      double precision:: finvl(nvar), finvr(nvar)
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,:)
c            btt(ilo:ihi)=Btaxis/REta(i,ilo:ihi)
            btt(ilo-1:ihi)=Btaxis/RC(i,ilo-1:ihi)
            call ConstructLRStatesEta(vd,vl,vr)
            do j=ilo,ihi,1
c
c$$$               finvl(1)=-Reta(i,j)*dZdXi(i,j)*(vl(j,4)+Btt(j)*vl(j,3))
c$$$     &              -Reta(i,j)*vl(j,1)*
c$$$     &              (-dZdXi(i,j)*vl(j,1)+dRdXi(i,j)*vl(j,2))
c$$$               finvl(2)=Reta(i,j)*dRdXi(i,j)*(vl(j,4)+Btt(j)*vl(j,3))
c$$$     &              -Reta(i,j)*vl(j,2)*
c$$$     &              (-dZdXi(i,j)*vl(j,1)+dRdXi(i,j)*vl(j,2))
c$$$               finvl(3)=0.D0
c$$$               finvl(4)=0.D0
c$$$c
c$$$               finvr(1)=-Reta(i,j)*dZdXi(i,j)*(vr(j,4)+Btt(j)*vr(j,3))
c$$$     &              -Reta(i,j)*vr(j,1)*
c$$$     &              (-dZdXi(i,j)*vr(j,1)+dRdXi(i,j)*vr(j,2))
c$$$               finvr(2)=Reta(i,j)*dRdXi(i,j)*(vr(j,4)+Btt(j)*vr(j,3))
c$$$     &              -Reta(i,j)*vr(j,2)*
c$$$     &              (-dZdXi(i,j)*vr(j,1)+dRdXi(i,j)*vr(j,2))
c$$$               finvr(3)=0.D0
c$$$               finvr(4)=0.D0
c     
               finvl(1)=-RC(i,j-1)*dZdXiC(i,j-1)*
     &              (vl(j,4)+Btt(j-1)*vl(j,3))
     &              -RC(i,j-1)*vl(j,1)*
     &              (-dZdXiC(i,j-1)*vl(j,1)+dRdXiC(i,j-1)*vl(j,2))
               finvl(2)=RC(i,j-1)*dRdXiC(i,j-1)*
     &              (vl(j,4)+Btt(j-1)*vl(j,3))
     &              -RC(i,j-1)*vl(j,2)*
     &              (-dZdXiC(i,j-1)*vl(j,1)+dRdXiC(i,j-1)*vl(j,2))
               finvl(3)=-Rc(i,j-1)*vl(j,3)*
     &              (-dZdXiC(i,j-1)*vl(j,1)+dRdXiC(i,j-1)*vl(j,2))
               finvl(4)=0.D0
c
               finvr(1)=-RC(i,j)*dZdXiC(i,j)*(vr(j,4)+Btt(j)*vr(j,3))
     &              -RC(i,j)*vr(j,1)*
     &              (-dZdXiC(i,j)*vr(j,1)+dRdXiC(i,j)*vr(j,2))
               finvr(2)=RC(i,j)*dRdXiC(i,j)*(vr(j,4)+Btt(j)*vr(j,3))
     &              -RC(i,j)*vr(j,2)*
     &              (-dZdXiC(i,j)*vr(j,1)+dRdXiC(i,j)*vr(j,2))
               finvr(3)=-RC(i,j)*vr(j,3)*
     &              (-dZdXiC(i,j)*vr(j,1)+dRdXiC(i,j)*vr(j,2))
               finvr(4)=0.D0
c
               finv(i,j,1:4)=0.5D0*(finvl(1:4)+finvr(1:4))
c$$$               write(61,*) i,j,btt(j),Rc(i,j),finv(i,j,1),
c$$$     &              finv(i,j,2),finvl(1),finvl(2), finvr(1),finvr(2)
c$$$               write(62,*) i,j,vl(j,1),vl(j,2),vl(j,3),vl(j,4),
c$$$     &              vr(j,1),vr(j,2),vr(j,3),vr(j,4)
c$$$               write(63,*) i,j-1,btt(j-1),Rc(i,j-1)
            enddo
         enddo
c     
      return
      end
c
c-----------------------------------------------------------------------