C -*- Mode: Fortran; -*-
c-----------------------------------------------------------------
c       Ravi Samtaney
c       Copyright 2004
c       Princeton Plasma Physics Laboratory
c       All Rights Reserved
c-----------------------------------------------------------------
c       $Log$
c-----------------------------------------------------------------
	subroutine SetBoundaryValuesRZ(Rn,Zn)
	use mesh_parms
	use mesh_common
	implicit none
	double precision:: Zn(ixlo:ixhi+1,iylo:iyhi+1)
	double precision:: Rn(ixlo:ixhi+1,iylo:iyhi+1)
	double precision:: Rma, Zma
	integer:: i,j,n
c       
c     Periodic BC in eta
c       If using Outer wall fixed
	j=0
	   do i=IXLO, IXHI+1,1
	      do n=0,nghost-1,1
		    Rn(i,j-n)=Rn(i,nylocal-n)
		    Zn(i,j-n)=Zn(i,nylocal-n)
	      enddo
	   enddo
	j=ny+1
	   do i=IXLO, IXHI,1
	      do n=0,nghost,1
		    Rn(i,j+n)=Rn(i,1+n)
		    Zn(i,j+n)=Zn(i,1+n)
	      enddo
	   enddo
c$$$c       If using magnetic axis fixed
c$$$	j=0
c$$$	   do i=IXLO, IXHI+1,1
c$$$	      do n=0,nghost-1,1
c$$$		    Rn(i,j-n)=Rn(i,nylocal-n)
c$$$		    Zn(i,j-n)=Zn(i,nylocal-n)
c$$$	      enddo
c$$$	   enddo
c$$$	j=ny+1
c$$$	   do i=IXLO, IXHI,1
c$$$		    Rn(i,1)=Rn(i,j)
c$$$		    Zn(i,1)=Zn(i,j)
c$$$	   enddo
c$$$	   do i=IXLO, IXHI,1
c$$$	      do n=0,nghost-1,1
c$$$		    Rn(i,j+n)=Rn(i,1+n)
c$$$		    Zn(i,j+n)=Zn(i,1+n)
c$$$	      enddo
c$$$	   enddo
c
c       Magnetic axis
c	   do j=1,nylocal+1,1
c	   i=1
c	   Rma=sum(Rn(i,1:nylocal+1))/(nylocal+1)
c	   Zma=sum(Zn(i,1:nylocal+1))/(nylocal+1)
c	   Rn(i,:)=Rma
c	   Zn(i,:)=Zma
	return
	end
c
c-----------------------------------------------------------------
	subroutine SetBoundaryValues(u)
	use mesh_parms
	use mesh_common
	use MagAxisToroidalField
	implicit none
	double precision:: u(IXLO:IXHI,IYLO:IYHI,nvar)
	integer i,j,l,n
	integer xbdry_reflect, ybdry_reflect, zbdry_reflect
c       
c       Boundary conditions are periodic in phi direction
c       For parallel, rely on MPI to make this happen.
c       For serial, it is done here.
c#ifndef PARALLEL
#ifdef ETAPERIODIC
	if(ny.eq.nylocal) then
	j=0
	   do i=IXLO, IXHI,1
	      do n=0,nghost-1,1
		 do l=1,nvar,1
		    u(i,j-n,l)=u(i,nylocal-n,l)
		 enddo
	      enddo
	   enddo
	j=ny+1
	   do i=IXLO, IXHI,1
	      do n=0,nghost-1,1
		 do l=1,nvar,1
		    u(i,j+n,l)=u(i,1+n,l)
		 enddo
	      enddo
	   enddo
	endif
#endif
	return
	end
c




c-----------------------------------------------------------------
	subroutine FluxBoundaryConditions(vx,finv,hinv)
	use mesh
	use mesh_common
	use BoundaryPressTemp
	use MagAxisToroidalField
	implicit none
	double precision:: vx(ixlo:ixhi,iylo:iyhi,nvar)
	double precision:: finv(1:nxlocal+1,1:nylocal,nvar)
	double precision:: hinv(1:nxlocal,1:nylocal+1,nvar)
c
c       
	double precision:: vl(nvar), vr(nvar), vm(nvar)
	double precision:: wl(nvar), wr(nvar)
	double precision:: fn(nvar), fxi(nvar)
	double precision:: costheta, sintheta
	double precision:: btt(1)
	integer:: i,j
c       
c       Ideal conductor boundary conditions
c       Left Boundary
	if(hasxlo.eq.1) then
c	   finv(1,:,3)=0.D0
	   finv(1,:,3:nvar)=0.D0
	   i=1
	   do j=1,nylocal,1
		 finv(i,j,1:2)=0.D0
	      enddo
	endif
c       Right Boundary
	if(hasxhi.eq.1) then
c	   finv(nxlocal+1,:,1)=0.D0
	   finv(nxlocal+1,:,3:nvar)=0.D0
	   i=nxlocal+1
	      do j=1,nylocal,1
		 vl(:)=vx(i-1,j,:)
		 btt(1)=Btaxis/RC(i-1,j)
		 fxi(1)=dZdeta(i,j)*Rxi(i,j)*(vl(4)+btt(1)*vl(3))
		 fxi(2)=-dRdeta(i,j)*Rxi(i,j)*(vl(4)+btt(1)*vl(3))
		 finv(i,j,1:2)=fxi(1:2)
	      enddo
	endif
c
	return
	end

c$$$c-----------------------------------------------------------------
c$$$	subroutine FluxBoundaryConditionsOld(vx,finv,hinv)
c$$$	use mesh
c$$$	use mesh_common
c$$$	use BoundaryPressTemp
c$$$	use MagAxisToroidalField
c$$$	implicit none
c$$$	double precision:: vx(ixlo:ixhi,iylo:iyhi,nvar)
c$$$	double precision:: finv(1:nxlocal+1,1:nylocal,nvar)
c$$$	double precision:: hinv(1:nxlocal,1:nylocal+1,nvar)
c$$$c
c$$$c       
c$$$	double precision:: vl(nvar), vr(nvar), vm(nvar)
c$$$	double precision:: wl(nvar), wr(nvar)
c$$$	double precision:: fn(nvar), fxi(nvar)
c$$$	double precision:: costheta, sintheta
c$$$	double precision:: btt(1)
c$$$	integer:: i,j
c$$$c       
c$$$c       Ideal conductor boundary conditions
c$$$c       Left Boundary
c$$$	if(hasxlo.eq.1) then
c$$$c	   finv(1,:,3)=0.D0
c$$$	   finv(1,:,3:nvar)=0.D0
c$$$	   i=1
c$$$	   do j=1,nylocal,1
c$$$#ifndef ETAPERIODIC
c$$$		 vr(:)=vx(i,j,:)
c$$$		 costheta=dZdEta(i,j)/dlXi(i,j)
c$$$		 sintheta=-dRdEta(i,j)/dlXi(i,j)
c$$$		 wr(1)=costheta*vr(1)+sintheta*vr(2)
c$$$		 wr(2)=-sintheta*vr(1)+costheta*vr(2)
c$$$		 wr(3:4)=vr(3:4)
c$$$		 wl(1)=0.D0
c$$$		 wl=wr
c$$$		 wl(1)=-wr(1)
c$$$		 btt(1)=Btaxis/Rxi(i,j)
c$$$		 call InviscidFluxRP(fn,vm,wl,wr,btt,1,1,1,1,nvar,1)
c$$$		 fxi(1)=fn(1)*costheta-
c$$$     &                          fn(2)*sintheta
c$$$		 fxi(2)=fn(1)*sintheta+
c$$$     &                          fn(2)*costheta
c$$$		 finv(i,j,1:2)=Rxi(i,j)*fxi(1:2)*dlXi(i,j)
c$$$#else
c$$$		 finv(i,j,1:2)=0.D0
c$$$#endif
c$$$	      enddo
c$$$	endif
c$$$c       Right Boundary
c$$$	if(hasxhi.eq.1) then
c$$$c	   finv(nxlocal+1,:,1)=0.D0
c$$$	   finv(nxlocal+1,:,3:nvar)=0.D0
c$$$	   i=nxlocal+1
c$$$	      do j=1,nylocal,1
c$$$		 vl(:)=vx(i-1,j,:)
c$$$		 costheta=dZdEta(i,j)/dlXi(i,j)
c$$$		 sintheta=-dRdEta(i,j)/dlXi(i,j)
c$$$		 wl(1)=costheta*vl(1)+sintheta*vl(2)
c$$$		 wl(1)=0.D0
c$$$		 wl(2)=-sintheta*vl(1)+costheta*vl(2)
c$$$		 wl(3:4)=vl(3:4)
c$$$		 wr=wl
c$$$		 wr(1)=-wl(1)
c$$$		 btt(1)=Btaxis/Rxi(i,j)
c$$$		 call InviscidFluxRP(fn,vm,wl,wr,btt,1,1,1,1,nvar,1)
c$$$		 fxi(1)=fn(1)*costheta-
c$$$     &                          fn(2)*sintheta
c$$$		 fxi(2)=fn(1)*sintheta+
c$$$     &                          fn(2)*costheta
c$$$		 finv(i,j,1:2)=Rxi(i,j)*fxi(1:2)*dlXi(i,j)
c$$$	      enddo
c$$$	endif
c
c       Bottom Boundary
c$$$	if(hasylo.eq.1) then
c$$$	   hinv(:,1,:,1)=0.D0
c$$$	   hinv(:,1,:,4:nvar)=0.D0
c$$$	   j=1
c$$$	      do i=1,nxlocal,1
c$$$		 vr(:)=vx(i,j,:)
c$$$		 costheta=-dZdXi(i,j)/dlEta(i,j)
c$$$		 sintheta=dRdXi(i,j)/dlEta(i,j)
c$$$		 wr(1)=vr(1)
c$$$		 wr(2)=costheta*vr(2)+sintheta*vr(3)
c$$$		 wr(3)=-sintheta*vr(2)+costheta*vr(3)
c$$$		 wr(4)=vr(4)
c$$$		 wr(5)=costheta*vr(5)+sintheta*vr(6)
c$$$	         wr(5)=0.D0
c$$$		 wr(6)=-sintheta*vr(5)+costheta*vr(6)
c$$$		 wr(7:8)=vr(7:8)
c$$$		 wl=wr
c$$$		 wl(2)=-wr(2)
c$$$		 wl(5)=-wr(5)
c$$$		 btt(1)=Btaxis/Reta(i,j)
c$$$		 call InviscidFluxRP(fn,vm,wl,wr,btt,1,1,1,1,nvar,2)
c$$$		 fxi(2)=fn(2)*costheta-
c$$$     &                          fn(3)*sintheta
c$$$		 fxi(3)=fn(2)*sintheta+
c$$$     &                          fn(3)*costheta
c$$$		 hinv(i,j,2:3)=Reta(i,j)*fxi(2:3)*dlEta(i,j)
c$$$	      enddo
c$$$	endif
c$$$c       Top Boundary
c$$$	if(hasyhi.eq.1) then
c$$$	   hinv(:,nylocal+1,:,1)=0.D0
c$$$	   hinv(:,nylocal+1,:,4:nvar)=0.D0
c$$$	   j=nylocal+1
c$$$	      do i=1,nxlocal,1
c$$$		 vl(:)=vx(i,j-1,:)
c$$$		 costheta=-dZdXi(i,j)/dlEta(i,j)
c$$$		 sintheta=dRdXi(i,j)/dlEta(i,j)
c$$$		 wl(1)=vl(1)
c$$$		 wl(2)=costheta*vl(2)+sintheta*vl(3)
c$$$		 wl(3)=-sintheta*vl(2)+costheta*vl(3)
c$$$		 wl(4)=vl(4)
c$$$		 wl(5)=costheta*vl(5)+sintheta*vl(6)
c$$$		 wl(6)=-sintheta*vl(5)+costheta*vl(6)
c$$$		 wl(7:8)=vl(7:8)
c$$$	         wl(5)=0.D0
c$$$		 wr=wl
c$$$		 wr(2)=-wl(2)
c$$$		 wr(5)=-wl(5)
c$$$		 btt(1)=Btaxis/Reta(i,j)
c$$$		 call InviscidFluxRP(fn,vm,wl,wr,btt,1,1,1,1,nvar,2)
c$$$		 fxi(2)=fn(2)*costheta-
c$$$     &                          fn(3)*sintheta
c$$$		 fxi(3)=fn(2)*sintheta+
c$$$     &                          fn(3)*costheta
c$$$		 hinv(i,j,2:3)=Reta(i,j)*fxi(2:3)*dlEta(i,j)
c$$$	      enddo
c$$$	endif
c$$$	return
c$$$	end