c
c
c =========================================================
      subroutine src2(maxmx,maxmy,meqn,mbc,mx,my,xlower,ylower,
     &                 dx,dy,q,maux,aux,t,dt)
c =========================================================
      implicit double precision(a-h,o-z)
      dimension    q(1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc, meqn)
      dimension  aux(1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc, maux)
c
c     # source terms for cylindrical symmetry in 2d Euler equations
c     # about y-axis (so x=radius)
c
c     # aux(i,j,1) stores x coordinate of cell center
c
      dimension qstar(4)
      common /param/  gamma,gamma1
c
c     # 2-stage Runge-Kutta method
c
      dt2    = dt/2.d0
      press  = 0.d0
      ndim   = 2

      do 10 i=1,mx
       do 10 j=1,my
         x        = aux(i,j,8)
         if (x .lt. 0) then
            write(6,'(A)') 'src2_axis : Radius (aux(i,j,8)) < 0'
            stop
         endif
	 rho      = q(i,j,1)
	 u        = q(i,j,2)/q(i,j,1)
	 v        = q(i,j,3)/q(i,j,1)
         e        = q(i,j,4)
 	 press    = gamma1*(e - 0.5d0*rho*(u**2 + v**2))

	 if (x.eq.0.d0) write(6,*) 'x = 0 in src2'
	 if (rho.eq.0.d0) write(6,*) 'rho = 0 in src2'

	 qstar(1) = q(i,j,1) - dt2*(ndim-1)/x * rho*u
	 qstar(2) = q(i,j,2) - dt2*(ndim-1)/x * rho*u**2
 	 qstar(3) = q(i,j,3) - dt2*(ndim-1)/x * rho*u*v
	 qstar(4) = q(i,j,4) - dt2*(ndim-1)/x * u*(e + press)
c
c        # second stage
c
	 rho      = qstar(1)
	 u        = qstar(2)/qstar(1)
	 v        = qstar(3)/qstar(1)
         e        = qstar(4)
 	 press    = gamma1*(e - 0.5d0*rho*(u**2 + v**2))
	 if (rho.eq.0.d0) write(6,*) 'rho = 0 in qstar'

	 q(i,j,1) = q(i,j,1) - dt*(ndim-1)/x * rho*u
	 q(i,j,2) = q(i,j,2) - dt*(ndim-1)/x * rho*u**2
 	 q(i,j,3) = q(i,j,3) - dt*(ndim-1)/x * rho*u*v
	 q(i,j,4) = q(i,j,4) - dt*(ndim-1)/x * u*(e + press)
   10    continue

      return
      end
