subroutine sod_initial

     	implicit none

!-------------------------------------------------------------------------------
! Include required header files. 
! Include files required for amr package
        use paramesh_dimensions
        use physicaldata
	use tree

! include files required for MPI library.
        include 'mpif.h'

! Pointer declarations
	include 'pointers.fh'


	integer	:: lblock, ierr
	logical	:: l_bc

!-------------------------------------------------------------------------------
! USER CODE BLOCK

	parameter(neq=nvar)					!<<< EDITED
	parameter(mx=nxb+2*nguard)				!<<< EDITED
	parameter(mxp=mx+1)					!<<< EDITED

	integer :: mx,mxp,nx, neq, nxp
        integer :: i

       	real, dimension(neq,mx) :: vg ! primitive variables

	common/minima/rmin,pmin,smin
	common/control/crn,fact,nriem,ibc

    	real    :: crn,fact,rmin,pmin,smin
	integer	:: nriem
        integer :: ibc



! Local Variables
	real	:: xpos,xjump
      	real    :: xmin,xmax, delx
        real    :: rr,rl,vr,vl, pr, pl


! Mesh Variables
        real, dimension(mx)    ::  rvol ! inverse control volume at i,j
	real, dimension(mx)    ::  xc
	real, dimension(mxp)   ::  xm   ! xm(i) is the cell interface position 
                                        ! to the left of the cell center
	real, dimension(mx)    ::  dxm
	real, dimension(mx)    ::  rdx
	real, dimension(mx)    ::  dxc  ! dx is the distance from xm(i) to 
                                        ! xm(i+1)
	real, dimension(mx)    ::  rdxc ! rdxc is the inverse of dxc



! Set advection algorithm control parameters
        crn   = 0.2			! courant number
        rmin  = 1.0e-5			! minimum allowed density
        pmin  = 1.0e-7			! minimum allowed pressure
        smin  = 1.0e-10
        fact  = 1.0
        nriem = 5			! number of iterations inside Riemann
                                        ! solver


	nx = nxb+2*nguard					!<<< EDITED
	nxp = nx+1


! END OF USER CODE BLOCK
!-------------------------------------------------------------------------------


! Loop over grid blocks 
	if(lnblocks.gt.0) then
	do lblock = 1,lnblocks


! use pointer to connect user variable to the amr package's solution array unk.
	uug => unk(:,:,lblock)

	xmax   = coord(1,1,lblock) + .5*bsize(1,lblock)
	xmin   = coord(2,1,lblock) - .5*bsize(1,lblock)
	delx   = (xmax - xmin)/float(nxb)


! Is this block at an external boundary?
	l_bc = .false.
	if(neigh(1,1,lblock).le.-20. .or.
     .			neigh(1,1,lblock).le.-20) l_bc = .true.


!-------------------------------------------------------------------------------
! USER CODE BLOCK
!
! This block of code is the users code for initializing the solution
! on an individual grid block.
!
!
! 1d Sod test problem
!
!     boundary condition factors
!
        ibc     = 2

!
!     create grid
!
      call ngrid ( xmin,xmax,delx,rdx,  
     .                   xc  ,xm  ,dxc ,dxm,  
     .                   rdxc,fact,rvol,
     .                   mx  ,mxp ,nx  , nxp   ) 
                    
!
! Sod problem initialization
!
	rr      =  0.125		! right side density
	rl      =  1.0			! left side density
	vr      =  0.0			! right side velocity
	vl      =  0.0			! left side velocity
	pr      =  0.1			! right side pressure
	pl      =  1.0			! left side pressure
        xjump   = .42			! location of diaphraghm
 

	do i=1,nx
       		vg(1,i) = rl
       		vg(2,i) = 0.0
       		vg(3,i) = pl
       		if (xc(i) .gt. xjump ) then
       			vg(1,i) = rr
       			vg(3,i) = pr
       		endif
	enddo



	call vgtouu(uug,vg,mx,nx,neq)

	if(l_bc) call bcset(uug,mx,nx,neq,rmin,pmin,ibc)

! END OF USER CODE BLOCK
!-------------------------------------------------------------------------------



	enddo
	endif

	call MPI_BARRIER(MPI_COMM_WORLD, ierr)



	return
	end subroutine sod_initial