c =========================================================
       subroutine qinit(maxmx,maxmy, meqn,mbc,mx,my, xlower,
     &      ylower, dx,dy, q,maux,aux)
c =========================================================
c
c     # Set initial conditions for q.

      implicit none

      integer maxmx, maxmy, meqn, mbc, mx, my, maux
      double precision xlower, ylower, dx, dy
      double precision q(1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc,meqn)
      double precision aux(1-mbc:maxmx+mbc,1-mbc:maxmy+mbc,maux)
      double precision xp,yp

      integer i, j, k, m
      double precision xlow, ylow, xc,yc, x0,y0, r0
      double precision p1, v1, rho1, z1, e1, rhou1, rhov1
      double precision p2, v2, rho2, z2, e2, rhou2, rhov2
      double precision press, rhou, rhov, rho,z, e
      double precision u2v21,u2v22, th, w
      double precision gamma, gamma1, fdisc


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

      e1 = p1/gamma1
      e2 = p2/gamma1


      do j = 1-mbc,my+mbc
         ylow = ylower + (j-1)*dy
         yc = ylower + (j-0.5d0)*dy
         do i = 1-mbc,mx+mbc
            xlow = xlower + (i-1)*dx
            xc = xlower + (i-0.5d0)*dx
            call cellave(xlow,ylow,dx,dy,w)
            rho = w*rho2 + (1-w)*rho1
            e = w*e2 + (1-w)*e1
            q(i,j,1) = rho
            q(i,j,2) = 0.d0
            q(i,j,3) = 0.d0
            q(i,j,4) = e
         enddo
      enddo

      return
      end


      double precision function fdisc(x,y)
      implicit none

      double precision x,y, x0, y0, r0, r, xp, yp

      common /cominit3/ x0, y0, r0

      call mapc2p(x,y,xp,yp)

      r = sqrt((xp-x0)**2 + (yp-y0)**2)
      fdisc = r - r0
      end
