subroutine projec(uf,vf,us,vs,px,py,nx,ny,mxy,piter)
c*******************************************************************
c
c     correct the velocity field 
c
c     (uf,vf)      <-     generated velocity field
c     (us,vs)      <-     velocity correction before project
c     grad cor     <-     gradient component of (us,vs)
c
c     local arrays:
c     cor,res,div:        1-d arrays used in the multigrid package
c                         to store the solution, residue and the rhs.
c     (adx,ady):          grad of cor.
c
c*******************************************************************
      integer nx,ny,i,j,k,kk,nxm,nym,mxy,mxym
      include 'gridsz'
      logical piter
      double precision uf(nx,ny),vf(nx,ny),
     +                 us(nx,ny),vs(nx,ny),
     +                 px(nx,ny),py(nx,ny)
      double precision cor(mxym),res(mxym),div(mxym),ro(mxym)
      double precision adx(nxm,nym),ady(nxm,nym)
      double precision hhinv,dt,roa
      common/timest/dt
      common/hfract/hhinv
      common/density/roa(0:nxm+1,0:nym+1)
      data cor/mxym*0.D0/
      save cor
      k(i,j) = (nx+3)*(j+1)+i+2

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

c***********************************************************
c     obtain the field to be projected
      do 5 j=1,ny
      do 5 i=1,nx
         us(i,j) = (us(i,j)-uf(i,j))/dt
         vs(i,j) = (vs(i,j)-vf(i,j))/dt
   5  continue
c
      call formdv(us,vs,div,nx,ny)
c
      call multi2(div,cor,ro,res,nx,ny,mxy)
c
c     generate the gradient of cor
      do 10 j=1,ny
      do 10 i=1,nx
         adx(i,j) = (cor(k(i,j))   - cor(k(i-1,j)) +
     +               cor(k(i,j-1)) - cor(k(i-1,j-1)) )*hhinv
         ady(i,j) = (cor(k(i,j))   - cor(k(i,j-1)) +
     +               cor(k(i-1,j)) - cor(k(i-1,j-1)) )*hhinv
 10   continue

c***************************************************
c     update the pressure
      do 30 j=1,ny
      do 30 i=1,nx
         px(i,j) = px(i,j) + adx(i,j)
         py(i,j) = py(i,j) + ady(i,j)
 30   continue
c
c     update the velocity
      if (.not.piter) then
      do 40 j=1,ny
      do 40 i=1,nx
         uf(i,j) = uf(i,j) + (us(i,j) - adx(i,j)/roa(i,j))*dt
         vf(i,j) = vf(i,j) + (vs(i,j) - ady(i,j)/roa(i,j))*dt
 40   continue
      end if

      end

      subroutine formdv(u,v,div,nx,ny)
c********************************************************************
c
c     (u,v)            <-      velocity field
c     div              ->      h^2 div (u,v)
c     assuming homogeneous bc.
c
c********************************************************************
      integer i,j,nx,ny
      double precision u(nx,ny),v(nx,ny),
     +                 div(-1:nx+1,-1:ny+1)
      double precision h,hinv,h2,sdiv,divm,dt
      common/timest/dt
      common/meshsz/h,hinv
      h2 = 0.5D0*h

      do 10 j=1,ny-1
      do 10 i=1,nx-1
         div(i,j) = h2*(
     +              u(i+1,j+1) - u(i,j+1) +
     +              u(i+1,j  ) - u(i,j  ) +
     +              v(i+1,j+1) - v(i+1,j) +
     +              v(i,  j+1) - v(i,  j) )
 10   continue
      do 20 j=1,ny-1
         div(0,j)  = h2*( 
     +               u(1,j+1) + u(1,j) +
     +               v(1,j+1) - v(1,j) )
         div(nx,j) = h2*(
     +               -u(nx,j+1) - u(nx,j) +
     +                v(nx,j+1) - v(nx,j) )
 20   continue
      do 30 i=1,nx-1
         div(i,0)  = h2*(
     +               u(i+1,1) - u(i,1) +
     +               v(i+1,1) + v(i,1) )
         div(i,ny) = h2*(
     +               u(i+1,ny) - u(i,ny) -
     +               v(i+1,ny) - v(i,ny))
 30   continue
      div(0,0)   =  h2*(u(1,1)   + v(1,1))
      div(0,ny)  =  h2*(u(1,ny)  - v(1,ny))
      div(nx,0)  = -h2*(u(nx,1)  - v(nx,1))
      div(nx,ny) = -h2*(u(nx,ny) + v(nx,ny))
c
      sdiv = 0.D0
      divm = 0.D0
      do 110 j=1,ny
      do 110 i=1,nx
         divm = dmax1(divm,dabs(div(i,j)))
         sdiv = sdiv + div(i,j)*div(i,j)
 110  continue
      sdiv = dsqrt(sdiv)*hinv
      divm = divm*hinv*hinv
      print *,'l2 of div',sdiv,' divm=',divm

      end