c#######################################################################
c#
c#                Cadez Differencing Coefficients
c#
c#######################################################################
      subroutine psiresid(Presid,
     .                    Psi,PsiD,KK,PsiTau,PsiRes,LapX,LapY,Volume,
     .                    ext,int,fil,UseTau)
      implicit none
c-----------------------------------------------------------------------
      integer ext(4),int(4),fil(4),UseTau
      real*8 Presid
      real*8 Psi(ext(1):ext(3),ext(2):ext(4))
      real*8 PsiD(ext(1):ext(3),ext(2):ext(4))
      real*8 KK(ext(1):ext(3),ext(2):ext(4))
      real*8 PsiTau(ext(1):ext(3),ext(2):ext(4))
      real*8 PsiRes(ext(1):ext(3),ext(2):ext(4))
      real*8 LapX(ext(1):ext(3),ext(2):ext(4))
      real*8 LapY(ext(1):ext(3),ext(2):ext(4))
      real*8 Volume(ext(1):ext(3),ext(2):ext(4))
c-----------------------------------------------------------------------
      integer i,j,imin,imax,jmin,jmax
      real*8 k2
      real*8 ZERO,F1o8,HALF,ONE,TWO,THREE,FOUR
      parameter(ZERO = 0.D0,F1o8 = 1.D0/8.D0,HALF = 0.5D0,ONE = 1.D0,
     .          TWO = 2.D0,THREE = 3.D0,FOUR = 4.D0)
c
c-----------------------------------------------------------------------
      imin = int(1)
      jmin = int(2)
      imax = int(3)
      jmax = int(4)
      if(imin.eq.fil(1)) imin = imin + 1
      if(jmin.eq.fil(2)) jmin = jmin + 1
      if(imax.eq.fil(3)) imax = imax - 1
      if(jmax.eq.fil(4)) jmax = jmax - 1
c
c-----------------------------------------------------------------------
c-  Set residal values for Psi & Lapse
c- NOT including Tau
c- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c- Inner Boundary
      if(jmin.ne.int(2)) then
        j = jmin-1
        if(imin.ne.int(1)) then
          i = imin-1
          k2 = F1o8*Volume(i,j)*KK(i,j)
          PsiRes(i,j) = LapX(i,j)  *(PsiD(i,j) - PsiD(i+1,j)) +
     .                  LapY(i,j)  *(PsiD(i,j) - PsiD(i,j+1)) +
     .                  LapX(i,j-1)*
     .                  LapY(i,j-1)*(PsiD(i,j) + ONE        ) -
     .                  k2*(Psi(i,j)**(-7))
        end if
        do i = imin,imax
          k2 = F1o8*Volume(i,j)*KK(i,j)
          PsiRes(i,j) = LapX(i,j)  *(PsiD(i,j) - PsiD(i+1,j)) +
     .                  LapX(i-1,j)*(PsiD(i,j) - PsiD(i-1,j)) +
     .                  LapY(i,j)  *(PsiD(i,j) - PsiD(i,j+1)) +
     .                  LapX(i,j-1)*
     .                  LapY(i,j-1)*(PsiD(i,j) + ONE        ) -
     .                  k2*(Psi(i,j)**(-7))
        end do
        if(imax.ne.int(3)) then
          i = imax+1
          k2 = F1o8*Volume(i,j)*KK(i,j)
          PsiRes(i,j) = LapX(i-1,j)*(PsiD(i,j) - PsiD(i-1,j)) +
     .                  LapY(i,j)  *(PsiD(i,j) - PsiD(i,j+1)) +
     .                  LapX(i,j-1)*
     .                  LapY(i,j-1)*(PsiD(i,j) + ONE        ) -
     .                  k2*(Psi(i,j)**(-7))
        end if
      end if
c- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c- Interior
      do j = jmin,jmax
        if(imin.ne.int(1)) then
          i = imin-1
          k2 = F1o8*Volume(i,j)*KK(i,j)
          PsiRes(i,j) = LapX(i,j)  *(PsiD(i,j) - PsiD(i+1,j)) +
     .                  LapY(i,j)  *(PsiD(i,j) - PsiD(i,j+1)) +
     .                  LapY(i,j-1)*(PsiD(i,j) - PsiD(i,j-1)) -
     .                  k2*(Psi(i,j)**(-7))
        end if
        do i = imin,imax
          k2 = F1o8*Volume(i,j)*KK(i,j)
          PsiRes(i,j) = LapX(i,j)  *(PsiD(i,j) - PsiD(i+1,j)) +
     .                  LapX(i-1,j)*(PsiD(i,j) - PsiD(i-1,j)) +
     .                  LapY(i,j)  *(PsiD(i,j) - PsiD(i,j+1)) +
     .                  LapY(i,j-1)*(PsiD(i,j) - PsiD(i,j-1)) -
     .                  k2*(Psi(i,j)**(-7))
        end do
        if(imax.ne.int(3)) then
          i = imax+1
          k2 = F1o8*Volume(i,j)*KK(i,j)
          PsiRes(i,j) = LapX(i-1,j)*(PsiD(i,j) - PsiD(i-1,j)) +
     .                  LapY(i,j)  *(PsiD(i,j) - PsiD(i,j+1)) +
     .                  LapY(i,j-1)*(PsiD(i,j) - PsiD(i,j-1)) -
     .                  k2*(Psi(i,j)**(-7))
        end if
      end do
c- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c- Outer Boundary
      if(jmax.ne.int(4)) then
        j = jmax+1
        if(imin.ne.int(1)) then
          i = imin-1
          k2 = F1o8*Volume(i,j)*KK(i,j)
          PsiRes(i,j) = LapX(i,j)  *(PsiD(i,j) - PsiD(i+1,j)) +
     .                  LapY(i,j)  *(PsiD(i,j)              ) +
     .                  LapY(i,j-1)*(PsiD(i,j) - PsiD(i,j-1)) -
     .                  k2*(Psi(i,j)**(-7))
        end if
        do i = imin,imax
          k2 = F1o8*Volume(i,j)*KK(i,j)
          PsiRes(i,j) = LapX(i,j)  *(PsiD(i,j) - PsiD(i+1,j)) +
     .                  LapX(i-1,j)*(PsiD(i,j) - PsiD(i-1,j)) +
     .                  LapY(i,j)  *(PsiD(i,j)              ) +
     .                  LapY(i,j-1)*(PsiD(i,j) - PsiD(i,j-1)) -
     .                  k2*(Psi(i,j)**(-7))
        end do
        if(imax.ne.int(3)) then
          i = imax+1
          k2 = F1o8*Volume(i,j)*KK(i,j)
          PsiRes(i,j) = LapX(i-1,j)*(PsiD(i,j) - PsiD(i-1,j)) +
     .                  LapY(i,j)  *(PsiD(i,j)              ) +
     .                  LapY(i,j-1)*(PsiD(i,j) - PsiD(i,j-1)) -
     .                  k2*(Psi(i,j)**(-7))
        end if
      end if
c
c-----------------------------------------------------------------------
c- Compute sum of square of residuals
c- Include Tau if required
c-----------------------------------------------------------------------
      if(UseTau.eq.1) then
		  do j = int(2),int(4)
		    do i = int(1),int(3)
            PsiRes(i,j) = PsiRes(i,j) - PsiTau(i,j)
            Presid = Presid + PsiRes(i,j)**2
c           Presid = Presid + ABS(PsiRes(i,j))
          end do
        end do
      else
		  do j = int(2),int(4)
		    do i = int(1),int(3)
            Presid = Presid + PsiRes(i,j)**2
c           Presid = Presid + ABS(PsiRes(i,j))
          end do
        end do
      end if
c
      return
      end
