c =========================================================
       subroutine qinit(maxmx,maxmy,maxmz,meqn,mbc,mx,my, mz,
     &      xlower,ylower,zlower,dx,dy,dz,q,maux,aux)
c =========================================================
c
c     # Set initial conditions for q.
c     # smooth density perturbation of magnitude rhoamp centered at x0
c     # with p, u chosen so it is a simple 3-wave before shock formation
c
c
      implicit none
      integer maxmx, maxmy, maxmz, meqn, mbc, mx, my, mz,maux
      double precision xlower, ylower, zlower, dx, dy, dz
      double precision q(1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc,
     &      1-mbc:maxmz+mbc,meqn)
      double precision aux(1-mbc:maxmx+mbc,1-mbc:maxmy+mbc,
     &      1-mbc:maxmz+mbc,maux)

      integer i,j,k,m,mcapa
      double precision e1, e2, xlow, ylow, zlow, w
      double precision p1, v1, rho1, p2, v2, rho2
      double precision xp, yp, zp
      double precision gamma, gamma1

      common /param/ gamma, gamma1
      common /cominit1/ p1, v1, rho1
      common /cominit2/ p2, v2, rho2

      e1 = p1/gamma1
      e2 = p2/gamma1

      do k = 1-mbc,mz+mbc
         zlow = zlower + (k-1)*dz
         do j = 1-mbc,my+mbc
            ylow = ylower + (j-1)*dy
            do i = 1-mbc,mx+mbc
               xlow = xlower + (i-1)*dx
               call cellave3(xlow,ylow,zlow,dx,dy,dz,w)
               q(i,j,k,1) = w*rho2 + (1-w)*rho1
               q(i,j,k,2) = 0.d0
               q(i,j,k,3) = 0.d0
               q(i,j,k,4) = 0.d0
               q(i,j,k,5) = w*e2 + (1-w)*e1
            enddo
         enddo
      enddo

      return
      end
