subroutine safePsiFromW(psi, 
     *     w, nX, nY, nMax,
     *     deltaX, deltaY, eps)

      implicit none

      real*8 psi, w, deltaX, deltaY, eps
      real*8 wTmp
      integer nX, nY, nMax,i,j

      real*8 diag, diagx, diagy
      real*8 psiLBound, psiHBound
 
      dimension wTmp(0:nMax+1,0:nMax+1)
      dimension psi(0:nMax+1,0:nMax+1)
      dimension w(0:nMax+1,0:nMax+1)

      call psiFromW(psi, 
     *     w, wTmp, nX, nY, nMax,
     *     deltaX, deltaY, eps)


      
      return
      end

      subroutine psiFromW(psi, 
     *     w, wTmp, nX, nY, nMax,
     *     deltaX, deltaY, eps)

      implicit none

      real*8 psi, w, deltaX, deltaY, eps
      real*8 wTmp
      integer nX, nY, nMax,i,j

      real*8 diag, diagx, diagy
      real*8 psiLBound, psiHBound
      real*8 getTopPsiBound
      real*8 getBottomPsiBound

      dimension wTmp(0:nMax+1,0:nMax+1)
      dimension psi(0:nMax+1,0:nMax+1)
      dimension w(0:nMax+1,0:nMax+1)

      diagx = -1.d0/(deltaX*deltaX)
      diagy = -1.d0/(deltaY*deltaY)

      diag = -2.d0*diagx-2.d0*diagy


c     This crud is just to let users change
c     the psi boundary conditions
c     Actually making use of it is strongly
c     discouranged, unless the user has carefully
c     thought through all the boundary issues
c     And even for those, changing the omega
c     boundary conditions makes more sense

      psiHBound = 
     *     getTopPsiBound()
      psiLBound = 
     *     getBottomPsiBound()


      do i=0,nX+1
         psi(i,0) = psiLBound
         psi(i,nY+1) = psiHBound
      enddo

c      write (*,*) "Before backup assign"

      do j=0,nY+1
         do i=0,nX+1
            wTmp(i,j) = w(i,j)
         enddo
      enddo

c      write (*,*) "Before recursion call"
      call conjgrad_2d_lat(diag, diagx, diagy,
     *     psi, w, nX, nY, nMax, eps, .FALSE.)

c      write (*,*) "After recursion call"

      do j=0,nY+1
         do i=0,nX+1
            w(i,j) = wTmp(i,j)
         enddo
      enddo

c      write (*,*) "After tmp re-assign"

      return
      end

c     ccccccccccccccccccccccccccccccccccccc
c     Always pair a test routine with the 
c     actual routine
c     
c     ccccccccccccccccccccccccccccccccccccc
      subroutine wFromPsi(psi, 
     *     w, nX, nY, nMax,
     *     deltaX, deltaY, eps)

      implicit none

      real*8 psi, w, deltaX, deltaY, eps
      real*8 multLat2d
      integer nX, nY, nMax, i, j

      real*8 diag, diagx, diagy

      dimension psi(0:nMax+1,0:nMax+1)
      dimension w(0:nMax+1,0:nMax+1)

      diagx = -1.d0/(deltaX*deltaX)
      diagy = -1.d0/(deltaY*deltaY)

      diag = -2.d0*diagx-2.d0*diagy

      do j=1,nY
         do i=1,nX
            
            w(i,j) = 
     *           multLat2d(diag,diagx,diagy,
     *           nX,nY, nMax, psi,i,j)

         enddo
      enddo

      return
      end