subroutine sweepy
C
C This subroutine performs sweeps in the Y direction,
C looping over i by k rows.   
C
C-----------------------------------------------------------------------
C
C GLOBALS
C
      include 'global.h'
      include 'zone.h'
      include 'sweep.h'
C
C-----------------------------------------------------------------------
C
C LOCALS
C
      integer i, j, k, n
      real paracoff(9, maxsweep)
C
C-----------------------------------------------------------------------
      sweep = 'y'
C
C Add two ghost zones on each side of 1D arrays
      nmin = 3
      nmax = jmax + 2
C
C Get 1D coordinates in order to set parabolic coefficients
      do 2 j = 1, jmax
         n = j + 2
         xa   (n) = zya(j)
         xa0  (n) = zya(j)
         dx   (n) = zdy(j)
         dx0  (n) = zdy(j)
 2    continue
C
C Set coordinates in ghost zones before calling paraset
      dx(nmin-1) = dx(nmin)
      dx(nmin-2) = dx(nmin+1)
      xa(nmin-1) = xa(nmin)   - dx(nmin-1)
      xa(nmin-2) = xa(nmin-1) - dx(nmin-2)
      dx(nmax+1) = dx(nmax)
      dx(nmax+2) = dx(nmax-1)
      xa(nmax+1) = xa(nmax)   + dx(nmax)
      xa0(nmax+1)= xa(nmax+1)
      xa(nmax+2) = xa(nmax+1) + dx(nmax+1)
      xa0(nmax+2)= xa(nmax+2)
C
C Set up coefficients for parabola subroutine.
      call paraset( paracoff, dx, xa, nmin, nmax, ngeomy)
C
C Now Loop over each column...
      radius = 1.
      do 1 k = 1, kmax
      do 1 i = 1, imax
C
C Put state variables into 1D arrays, padding with 2 ghost zones
      do 10 j = 1, jmax
         n = j + 2
         rho (n) = zrho(i,j,k)
         p   (n) = zp  (i,j,k)
         u   (n) = zuy (i,j,k)
         v   (n) = zuz (i,j,k)
         w   (n) = zux (i,j,k)
         xa  (n) = zya (j)
         dx  (n) = zdy (j)
         enrg(n) = p(n)/(rho(n)*gamm) 
     &           + 0.5*(u(n)**2+v(n)**2+w(n)**2)
 10   continue
         xa(nmax+1) = xa(nmax) + dx(nmax)
C
C Calculate the volume element appropriate to the desired geometry
C
      if(ngeomy .eq. 0) then
         do 20 n = nmin, nmax+1
            dvol (n) = dx(n)
            dvol0(n) = dvol(n)
 20      continue
      else if(ngeomy .eq. 1) then
         do 21 n = nmin, nmax+1
            dvol (n) = dx(n) * (xa(n) + .5*dx(n))
            dvol0(n) = dvol(n)
 21      continue
      else if(ngeomy .eq. 3) then
         radius = zxa(i) + 0.5*zdx(i)
         do 23 n = nmin, nmax+1
            dvol (n) = dx(n) * radius
            dvol0(n) = dvol(n)
 23      continue
      else if(ngeomy .eq. 4) then
         radius = zxa(i) + 0.5*zdx(i)
         do 24 n = nmin, nmax+1
            dvol (n) = (cos(xa(n)) - cos(xa(n+1))) * radius
            dvol0(n) = dvol(n)
 24      continue
      endif
C
C First call sweepbc to impose boundary conditions on the 1D arrays,
C  then call ppm.f to evolve the flow in the 1D sweep,
C   then call sweepbc to reset boundary values for remap,
C    then call remap.f to remap evolved quantities back to Eulerian grid.
C
      call sweepbc( nlefty, nrighty )
      call ppm    ( nlefty, nrighty, ngeomy, paracoff )
      call sweepbc( nlefty, nrighty )
      call remap  ( nlefty, nrighty, ngeomy )
C
C Put updated values into 2D arrays, dropping ghost zones
C
      do 100 j = 1, jmax
         n = j + 2
         zrho(i,j,k) = rho(n)
         zp  (i,j,k) = p  (n)
         zuy (i,j,k) = u  (n)
         zuz (i,j,k) = v  (n)
         zux (i,j,k) = w  (n)
 100  continue
C
 1    continue
C
      return
      end