c#######################################################################
c#
c#                Cadez Differencing Coefficients
c#
c#######################################################################
      subroutine psiresid(Presid,
     .                    Psi,PsiD,KK,PsiTau,PsiRes,
     .                    LapW,LapX,LapY,Volume,
     .                    ext,int,fil,UseTau)
      implicit none
c-----------------------------------------------------------------------
      integer ext(6),int(6),fil(6),UseTau
      real*8 Presid
      real*8 Psi(ext(1):ext(4),ext(2):ext(5),ext(3):ext(6))
      real*8 PsiD(ext(1):ext(4),ext(2):ext(5),ext(3):ext(6))
      real*8 KK(ext(1):ext(4),ext(2):ext(5),ext(3):ext(6))
      real*8 PsiTau(ext(1):ext(4),ext(2):ext(5),ext(3):ext(6))
      real*8 PsiRes(ext(1):ext(4),ext(2):ext(5),ext(3):ext(6))
      real*8 LapW(ext(2):ext(5),ext(3):ext(6))
      real*8 LapX(ext(2):ext(5),ext(3):ext(6))
      real*8 LapY(ext(2):ext(5),ext(3):ext(6))
      real*8 Volume(ext(2):ext(5),ext(3):ext(6))
c-----------------------------------------------------------------------
      integer k,i,j,kmin,imin,imax,kmax,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-----------------------------------------------------------------------
      kmin = int(1)
      imin = int(2)
      jmin = int(3)
      kmax = int(4)
      imax = int(5)
      jmax = int(6)
      if(imin.eq.fil(2)) imin = imin + 1
      if(jmin.eq.fil(3)) jmin = jmin + 1
      if(imax.eq.fil(5)) imax = imax - 1
      if(jmax.eq.fil(6)) jmax = jmax - 1
c
c-----------------------------------------------------------------------
c-  Set residal values for Psi & Lapse
c- NOT including Tau
c- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c- Inner Boundary
      if(jmin.ne.int(3)) then
        j = jmin-1
        if(imin.ne.int(2)) then
          i = imin-1
c     write(*,*) kmin-1,i,j,PsiD(kmin-1,i,j),Psi(kmin-1,i,j)
          do k = kmin,kmax
            k2 = F1o8*Volume(i,j)*KK(k,i,j)
            PsiRes(k,i,j) = LapW(i,j)  *(TWO*
     .                                   PsiD(k,i,j) - PsiD(k+1,i,j)
     .                                               - PsiD(k-1,i,j)) +
     .                      LapX(i,j)  *(PsiD(k,i,j) - PsiD(k,i+1,j)) +
     .                      LapY(i,j)  *(PsiD(k,i,j) - PsiD(k,i,j+1)) +
     .                      LapX(i,j-1)*
     .                      LapY(i,j-1)*(PsiD(k,i,j) + ONE          ) -
     .                      k2*(Psi(k,i,j)**(-7))
c     write(*,*) k,i,j,PsiD(k,i,j),Psi(k,i,j)
          end do
c     write(*,*) kmax+1,i,j,PsiD(kmax+1,i,j),Psi(kmax+1,i,j)
        end if
        do i = imin,imax
c     write(*,*) kmin-1,i,j,PsiD(kmin-1,i,j),Psi(kmin-1,i,j)
          do k = kmin,kmax
            k2 = F1o8*Volume(i,j)*KK(k,i,j)
            PsiRes(k,i,j) = LapW(i,j)  *(TWO*
     .                                   PsiD(k,i,j) - PsiD(k+1,i,j)
     .                                               - PsiD(k-1,i,j)) +
     .                      LapX(i,j)  *(PsiD(k,i,j) - PsiD(k,i+1,j)) +
     .                      LapX(i-1,j)*(PsiD(k,i,j) - PsiD(k,i-1,j)) +
     .                      LapY(i,j)  *(PsiD(k,i,j) - PsiD(k,i,j+1)) +
     .                      LapX(i,j-1)*
     .                      LapY(i,j-1)*(PsiD(k,i,j) + ONE          ) -
     .                      k2*(Psi(k,i,j)**(-7))
c     write(*,*) k,i,j,PsiD(k,i,j),Psi(k,i,j)
          end do
c     write(*,*) kmax+1,i,j,PsiD(kmax+1,i,j),Psi(kmax+1,i,j)
        end do
        if(imax.ne.int(5)) then
          i = imax+1
c     write(*,*) kmin-1,i,j,PsiD(kmin-1,i,j),Psi(kmin-1,i,j)
          do k = kmin,kmax
            k2 = F1o8*Volume(i,j)*KK(k,i,j)
            PsiRes(k,i,j) = LapW(i,j)  *(TWO*
     .                                   PsiD(k,i,j) - PsiD(k+1,i,j)
     .                                               - PsiD(k-1,i,j)) +
     .                      LapX(i-1,j)*(PsiD(k,i,j) - PsiD(k,i-1,j)) +
     .                      LapY(i,j)  *(PsiD(k,i,j) - PsiD(k,i,j+1)) +
     .                      LapX(i,j-1)*
     .                      LapY(i,j-1)*(PsiD(k,i,j) + ONE          ) -
     .                      k2*(Psi(k,i,j)**(-7))
c     write(*,*) k,i,j,PsiD(k,i,j),Psi(k,i,j)
          end do
c     write(*,*) kmax+1,i,j,PsiD(kmax+1,i,j),Psi(kmax+1,i,j)
        end if
      end if
c- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c- Interior
      do j = jmin,jmax
        if(imin.ne.int(2)) then
          i = imin-1
c     write(*,*) kmin-1,i,j,PsiD(kmin-1,i,j),Psi(kmin-1,i,j)
          do k = kmin,kmax
            k2 = F1o8*Volume(i,j)*KK(k,i,j)
            PsiRes(k,i,j) = LapW(i,j)  *(TWO*
     .                                   PsiD(k,i,j) - PsiD(k+1,i,j)
     .                                               - PsiD(k-1,i,j)) +
     .                      LapX(i,j)  *(PsiD(k,i,j) - PsiD(k,i+1,j)) +
     .                      LapY(i,j)  *(PsiD(k,i,j) - PsiD(k,i,j+1)) +
     .                      LapY(i,j-1)*(PsiD(k,i,j) - PsiD(k,i,j-1)) -
     .                      k2*(Psi(k,i,j)**(-7))
c     write(*,*) k,i,j,PsiD(k,i,j),Psi(k,i,j)
          end do
c     write(*,*) kmax+1,i,j,PsiD(kmax+1,i,j),Psi(kmax+1,i,j)
        end if
        do i = imin,imax
c     write(*,*) kmin-1,i,j,PsiD(kmin-1,i,j),Psi(kmin-1,i,j)
          do k = kmin,kmax
            k2 = F1o8*Volume(i,j)*KK(k,i,j)
            PsiRes(k,i,j) = LapW(i,j)  *(TWO*
     .                                   PsiD(k,i,j) - PsiD(k+1,i,j)
     .                                               - PsiD(k-1,i,j)) +
     .                      LapX(i,j)  *(PsiD(k,i,j) - PsiD(k,i+1,j)) +
     .                      LapX(i-1,j)*(PsiD(k,i,j) - PsiD(k,i-1,j)) +
     .                      LapY(i,j)  *(PsiD(k,i,j) - PsiD(k,i,j+1)) +
     .                      LapY(i,j-1)*(PsiD(k,i,j) - PsiD(k,i,j-1)) -
     .                      k2*(Psi(k,i,j)**(-7))
c     write(*,*) k,i,j,PsiD(k,i,j),Psi(k,i,j)
          end do
c     write(*,*) kmax+1,i,j,PsiD(kmax+1,i,j),Psi(kmax+1,i,j)
        end do
        if(imax.ne.int(5)) then
          i = imax+1
c     write(*,*) kmin-1,i,j,PsiD(kmin-1,i,j),Psi(kmin-1,i,j)
          do k = kmin,kmax
            k2 = F1o8*Volume(i,j)*KK(k,i,j)
            PsiRes(k,i,j) = LapW(i,j)  *(TWO*
     .                                   PsiD(k,i,j) - PsiD(k+1,i,j)
     .                                               - PsiD(k-1,i,j)) +
     .                      LapX(i-1,j)*(PsiD(k,i,j) - PsiD(k,i-1,j)) +
     .                      LapY(i,j)  *(PsiD(k,i,j) - PsiD(k,i,j+1)) +
     .                      LapY(i,j-1)*(PsiD(k,i,j) - PsiD(k,i,j-1)) -
     .                      k2*(Psi(k,i,j)**(-7))
c     write(*,*) k,i,j,PsiD(k,i,j),Psi(k,i,j)
          end do
c     write(*,*) kmax+1,i,j,PsiD(kmax+1,i,j),Psi(kmax+1,i,j)
        end if
      end do
c- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c- Outer Boundary
      if(jmax.ne.int(6)) then
        j = jmax+1
        if(imin.ne.int(2)) then
          i = imin-1
c     write(*,*) kmin-1,i,j,PsiD(kmin-1,i,j),Psi(kmin-1,i,j)
          do k = kmin,kmax
            k2 = F1o8*Volume(i,j)*KK(k,i,j)
            PsiRes(k,i,j) = LapW(i,j)  *(TWO*
     .                                   PsiD(k,i,j) - PsiD(k+1,i,j)
     .                                               - PsiD(k-1,i,j)) +
     .                      LapX(i,j)  *(PsiD(k,i,j) - PsiD(k,i+1,j)) +
     .                      LapY(i,j)  *(PsiD(k,i,j)                ) +
     .                      LapY(i,j-1)*(PsiD(k,i,j) - PsiD(k,i,j-1)) -
     .                      k2*(Psi(k,i,j)**(-7))
c     write(*,*) k,i,j,PsiD(k,i,j),Psi(k,i,j)
          end do
c     write(*,*) kmax+1,i,j,PsiD(kmax+1,i,j),Psi(kmax+1,i,j)
        end if
        do i = imin,imax
c     write(*,*) kmin-1,i,j,PsiD(kmin-1,i,j),Psi(kmin-1,i,j)
          do k = kmin,kmax
            k2 = F1o8*Volume(i,j)*KK(k,i,j)
            PsiRes(k,i,j) = LapW(i,j)  *(TWO*
     .                                   PsiD(k,i,j) - PsiD(k+1,i,j)
     .                                               - PsiD(k-1,i,j)) +
     .                      LapX(i,j)  *(PsiD(k,i,j) - PsiD(k,i+1,j)) +
     .                      LapX(i-1,j)*(PsiD(k,i,j) - PsiD(k,i-1,j)) +
     .                      LapY(i,j)  *(PsiD(k,i,j)                ) +
     .                      LapY(i,j-1)*(PsiD(k,i,j) - PsiD(k,i,j-1)) -
     .                      k2*(Psi(k,i,j)**(-7))
c     write(*,*) k,i,j,PsiD(k,i,j),Psi(k,i,j)
          end do
c     write(*,*) kmax+1,i,j,PsiD(kmax+1,i,j),Psi(kmax+1,i,j)
        end do
        if(imax.ne.int(5)) then
          i = imax+1
c     write(*,*) kmin-1,i,j,PsiD(kmin-1,i,j),Psi(kmin-1,i,j)
          do k = kmin,kmax
            k2 = F1o8*Volume(i,j)*KK(k,i,j)
            PsiRes(k,i,j) = LapW(i,j)  *(TWO*
     .                                   PsiD(k,i,j) - PsiD(k+1,i,j)
     .                                               - PsiD(k-1,i,j)) +
     .                      LapX(i-1,j)*(PsiD(k,i,j) - PsiD(k,i-1,j)) +
     .                      LapY(i,j)  *(PsiD(k,i,j)                ) +
     .                      LapY(i,j-1)*(PsiD(k,i,j) - PsiD(k,i,j-1)) -
     .                      k2*(Psi(k,i,j)**(-7))
c     write(*,*) k,i,j,PsiD(k,i,j),Psi(k,i,j)
          end do
c     write(*,*) kmax+1,i,j,PsiD(kmax+1,i,j),Psi(kmax+1,i,j)
        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(3),int(6)
          do i = int(2),int(5)
            do k = int(1),int(4)
              PsiRes(k,i,j) = PsiRes(k,i,j) - PsiTau(k,i,j)
              Presid = Presid + PsiRes(k,i,j)**2
c             Presid = Presid + ABS(PsiRes(k,i,j))
            end do
          end do
        end do
      else
        do j = int(3),int(6)
          do i = int(2),int(5)
            do k = int(1),int(4)
              Presid = Presid + PsiRes(k,i,j)**2
c             Presid = Presid + ABS(PsiRes(k,i,j))
            end do
          end do
        end do
      end if
c
      return
      end
