subroutine macpro(ux,vx,uy,vy,nx,ny,mxy)
**********************************************************
*     correcting the edge velocities by subtracting the
*     pressure gradient from the mac projection.
*
*         ---(uy,vy)----
*                       |
*                       |
*             (i,j)  (ux,vx)
*                       |
*                       |
*         --------------
***********************************************************
      integer i,j,nx,ny,k,nxm,nym,mxym,mxy
      include 'gridsz'
      double precision ux(0:nx,0:ny),vx(0:nx,0:ny),
     +                 uy(0:nx,0:ny),vy(0:nx,0:ny)
      double precision adx(0:nxm,nym),ady(nxm,0:nym),
     +                 cor(mxym),res(mxym),div(mxym),ro(mxym)
      double precision h,hinv,roa,rr,rt
      common/meshsz/h,hinv
      common/density/roa(0:nxm+1,0:nym+1)
      data cor/mxym*0.D0/
      save cor

      k(i,j) = (nx+2)*j+i+1

      do 2 i=1,mxym
         res(i) = 0.D0
         div(i) = 0.D0
 2    continue
      do 3 j=0,ny+1
      do 3 i=0,nx+1
         ro(k(i,j)) = 1.D0/roa(i,j)
 3    continue

c     obtain the divergence for the edge velocity field (MAC projection),
c     multiplied by h^2
      call fmdved(ux,vx,uy,vy,div,nx,ny)

      call multi1(div,cor,ro,res,nx,ny,mxy)

      do 10 j=1,ny
      do 10 i=1,nx-1
         rr = 0.5D0*(ro(k(i,j))+ro(k(i+1,j)))
         adx(i,j) = hinv*(cor(k(i+1,j))-cor(k(i,j)))*rr
 10   continue
      do 12 j=1,ny
         adx(0,j) = 0.D0
         adx(nx,j) = 0.D0
 12   continue

      do 20 j=1,ny-1
      do 20 i=1,nx
         rt = 0.5D0*(ro(k(i,j))+ro(k(i,j+1)))
         ady(i,j) = hinv*(cor(k(i,j+1))-cor(k(i,j)))*rt
 20   continue
      do 22 i=1,nx
         ady(i,0) = 0.D0
         ady(i,ny) = 0.D0
 22   continue

      do 30 j=1,ny
      do 30 i=1,nx-1
         ux(i,j) = ux(i,j) - adx(i,j)
 30   continue
      do 40 j=1,ny-1
      do 40 i=1,nx
         vy(i,j) = vy(i,j) - ady(i,j)
 40   continue

      do 50 j=1,ny-1
      do 50 i=1,nx
         uy(i,j) = uy(i,j) - 0.25D0*(
     +             adx(i,j) + adx(i-1,j) +
     +             adx(i,j+1) + adx(i-1,j+1) )
 50   continue
      do 60 j=1,ny
      do 60 i=1,nx-1
         vx(i,j) = vx(i,j) - 0.25D0*(
     +             ady(i,j) + ady(i+1,j) +
     +             ady(i,j-1) + ady(i+1,j-1) )
 60   continue

      end

      subroutine fmdved(ux,vx,uy,vy,div,nx,ny)
      integer i,j,nx,ny
      logical prnorm
      double precision ux(0:nx,0:ny),vx(0:nx,0:ny),
     +                 uy(0:nx,0:ny),vy(0:nx,0:ny),
     +                 div(0:nx+1,0:ny+1)
      double precision divm,sdiv,h,hinv
      common/meshsz/h,hinv
      parameter (prnorm=.false.)
      do 5 j=1,ny
      do 5 i=1,nx
         div(i,j) = h*( ux(i,j) - ux(i-1,j) + 
     +                  vy(i,j) - vy(i,j-1) )
  5   continue

C     if (prnorm .eq. .false.) return
      if (.not. prnorm) return

      divm = 0.D0
      sdiv = 0.D0
      do 10 j=1,ny
      do 10 i=1,nx
         divm = dmax1(divm,dabs(div(i,j)))
         sdiv = sdiv + div(i,j)**2
 10   continue
      sdiv = dsqrt(sdiv)*hinv
      divm = divm*hinv*hinv
      print *,'mac l2=',sdiv,' l_inf=',divm
      end