c
c
c
c     =====================================================
       subroutine qinit(maxmx,maxmy,meqn,mbc,mx,my,xlower,ylower,
     &                   dx,dy,q,maux,aux)
c     =====================================================
c
c     # Set initial conditions for q.
c
      implicit double precision (a-h,o-z)
      dimension q(1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc, meqn)
      common /cqinit/ x1,x2,y1,y2
c
      xm = (x1 + x2) / 2.d0
      ym = (y1 + y2) / 2.d0
      xscale = datan(1.d0)*8.d0 / (x2-x1)
      yscale = datan(1.d0)*8.d0 / (y2-y1)

      do 20 i=1,mx
         xc = xlower + (i-0.5d0)*dx
         do 20 j=1,my
            yc = ylower + (j-0.5d0)*dy
            call mapc2p(xc,yc,xp,yp)
	    if (xp.gt.x1 .and. xp.lt.x2 .and. 
     &	        yp.gt.y1 .and. yp.lt.y2) then
	        
                q(i,j,1) = (1.d0 + dcos((xp-xm)*xscale)) *
     &                     (1.d0 + dcos((yp-ym)*yscale))
              else
                q(i,j,1) = 0.d0 
              endif
  20        continue
      return
      end