subroutine vortic(u,v,vor,nx,ny)
c********************************************************************
c
c     obtain the vorticity of the velocity field.
c
c     (u,v)          ->      given velocity field
c     vor              ->      vorticity field.
c
c********************************************************************
      integer i,j,nx,ny
      double precision u(nx,ny),v(nx,ny),vor(nx,ny)
      double precision hhinv,thinv
      common/hfract/hhinv
      thinv = 2.D0/3.D0*hhinv
      do 10 j=2,ny-1
      do 10 i=2,nx-1
         vor(i,j) = hhinv*( v(i+1,j) - v(i-1,j) -
     +                      u(i,j+1) + u(i,j-1))
 10   continue
      do 20 j=2,ny-1
         vor(1,j)  =  thinv*( v(2,j) + 3.D0*v(1,j) ) +
     +                hhinv*( u(1,j-1) - u(1,j+1) )
         vor(nx,j) = -thinv*( v(nx-1,j) + 3.D0*v(nx,j) ) +
     +                hhinv*( u(nx,j-1) - u(nx,j+1) )
 20   continue
      do 30 i=2,nx-1
         vor(i,1)  = hhinv*( v(i+1,1) - v(i-1,1) ) -
     +               thinv*( u(i,2) + 3.D0*u(i,1) )
         vor(i,ny) = hhinv*( v(i+1,ny) - v(i-1,ny) ) +
     +               thinv*( u(i,ny-1) + 3.D0*u(i,ny) )
 30   continue
      vor(1,1)   = thinv*( v(2,1) + 3.D0*v(1,1) -
     +                     u(1,2) - 3.D0*u(1,1) )
      vor(1,ny)  = thinv*( v(2,ny) + 3.D0*v(1,ny) +
     +                     u(1,ny-1) + 3.D0*u(1,ny) )
      vor(nx,1)  = thinv*( -v(nx-1,1) - 3.D0*v(nx,1) -
     +                      u(nx,2) - 3.D0*u(nx,1) )
      vor(nx,ny) = thinv*( -v(nx-1,ny) - 3.D0*v(nx,ny) +
     +                      u(nx,ny-1) + 3.D0*u(nx,ny) )
c
      write(24,*) vor
      end