c     This routine does not fill the Y boundaries
c     An external call must take care of that

      subroutine uFromPsi(u, psi, nX, nY, 
     *     nMax, deltaX, deltaY, eps)

      implicit none

      real*8 psi, u, deltaX, deltaY, eps

      real*8 oneOver2dX
      real*8 oneOver2dY

      real*8 dPsiDx, dPsiDy

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

      dimension psi(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     U = curl(psi e_z)
c     U[x] = d Psi / dy
c     U[y] = -d Psi / dx

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

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

            dPsiDx =
     *           (psi(iPlus1,j) -
     *           psi(iMinus1,j))*
     *           oneOver2dX

            dPsiDy =
     *           (psi(i,j+1) -
     *           psi(i,j-1))*
     *           oneOver2dY

            u(1,i,j) = dPsiDy
            u(2,i,j) = -dPsiDx

         enddo
      enddo


      return
      end