c     W is such that w e_z = curl u
c     This just calculates the z
c     component (actually the only component)
c     of curl u
c
c     SHould only be called once at 
c     the beginning of the simulation
c     
c     The similar uFromPsi function
c     should be called every time output is
c     needed
c     
c     Assumes a uniform grid
c     deltaX and deltaY are
c     the grid spacing
c     
      subroutine wFromU(w, u, nX, nY, 
     *     nMax, deltaX, deltaY, eps)

      implicit none

      real*8 w, u, deltaX, deltaY, eps
      real*8 dvdx, dudy
      real*8 oneOver2dX, oneOVer2dY

      real*8 topBound, bottomBound
      real*8 getTopWBound
      real*8 getBottomWBound

      integer nX, nY, nMax, i, j
      integer iPlus1, iMinus1

      dimension w(0:nMax+1, 0:nMax+1)
      dimension u(2, 0:nMax+1, 0:nMax+1)

      oneOver2dX = 2.d0*deltaX
      oneOver2dX = 1.d0/oneOver2dX

      oneOver2dY = 2.d0*deltaY
      oneOver2dY = 1.d0/oneOver2dY

c     Just to give user option of changing
c     boundaries conditions on omega
c
c     Not reccomended unless user has
c     really thought things through
c
      topBound = getTopWBound()
      bottomBound = getBottomWBound()


c     w e_z = curl u
c     curl u = <0 , 0 , dv/dx - du/dy >

      do j=1,nY
         do i=1,nX

            iPlus1 = mod(i, nX)+1
            iMinus1 = mod(nX+i-2, nX)+1

            dvdx = (u(2,iPlus1,j)
     *           -u(2,iMinus1,j))
     *           *oneOVer2dX


            dudy = (u(1,i,j+1)
     *           -u(1,i,j-1))
     *           *oneOVer2dY

            w(i,j) = dvdx-dudy

            
c            write (*,*) i,j,dvdx,dudy,w(i,j), 
c     *           iMinus1, iPlus1
c            pause "i,j,dvdx,dudy,w(i,j)"
            
         end do
      end do

      do i=1,nX
         w(i,0) = topBound
         w(i,nY+1) = bottomBound
      end do

      return
      end