c#######################################################################
c#
c#                Cadez Smoothers
c#
c#######################################################################
c
c#######################################################################
c# Psi Red/Black Smoother
c#######################################################################
      subroutine psismoothrb(Psi,PsiD,KK,
     .                       PsiTau,LapX,LapY,Volume,
     .                       ext,int,fil,redblack)
      implicit none
c-----------------------------------------------------------------------
      integer ext(4),int(4),fil(4),redblack
      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 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,is
      real*8 numpsi, denpsi, k2
      integer red,black
      parameter(red = 0,black = 1)
      real*8 F1o8,ONE,TWO,THREE,FOUR,SIX,SEVEN
      parameter(F1o8 = 1.D0/8.D0,ONE = 1.D0,TWO = 2.D0,
     .          THREE = 3.D0,FOUR = 4.D0,SIX = 6.D0,SEVEN = 7.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- Perform Red/Black checkerboard smoothing for Psi
c- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c- Inner Boundary
      if(jmin.ne.int(2)) then
        j = jmin-1
        is = ABS(mod(j+redblack,2))
        do i = int(1)+is,int(3),2
          k2 = F1o8*Volume(i,j)*(Psi(i,j)**(-8))*KK(i,j)
          if((i.eq.int(1)).and.(imin.ne.int(1))) then
            numpsi = LapX(i,j)  *(PsiD(i+1,j) - PsiD(i,j)) +
     .               LapY(i,j)  *(PsiD(i,j+1) - PsiD(i,j)) -
     .               LapX(i,j-1)*
     .               LapY(i,j-1)*(ONE         + PsiD(i,j)) +
     .               Psi(i,j)*k2
            denpsi = ONE/(LapX(i,j) +
     .                    LapY(i,j) + LapX(i,j-1)*LapY(i,j-1) +
     .                    SEVEN*k2)
          else if((i.eq.int(3)).and.(imax.ne.int(3))) then
            numpsi = LapX(i-1,j)*(PsiD(i-1,j) - PsiD(i,j)) +
     .               LapY(i,j)  *(PsiD(i,j+1) - PsiD(i,j)) -
     .               LapX(i,j-1)*
     .               LapY(i,j-1)*(ONE         + PsiD(i,j)) +
     .               Psi(i,j)*k2
            denpsi = ONE/(              LapX(i-1,j) +
     .                    LapY(i,j)   + LapX(i,j-1)*LapY(i,j-1) +
     .                    SEVEN*k2      )
          else
            numpsi = LapX(i,j)  *(PsiD(i+1,j) - PsiD(i,j)) +
     .               LapX(i-1,j)*(PsiD(i-1,j) - PsiD(i,j)) +
     .               LapY(i,j)  *(PsiD(i,j+1) - PsiD(i,j)) -
     .               LapX(i,j-1)*
     .               LapY(i,j-1)*(ONE         + PsiD(i,j)) +
     .               Psi(i,j)*k2
            denpsi = ONE/(LapX(i,j) + LapX(i-1,j) +
     .                    LapY(i,j) + LapX(i,j-1)*LapY(i,j-1) +
     .                    SEVEN*k2)
          end if
          PsiD(i,j) = PsiD(i,j) + (numpsi + PsiTau(i,j))*denpsi
          Psi(i,j) = PsiD(i,j) + ONE
        end do
      end if
c- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c- Interior
      do j = jmin,jmax
        is = ABS(mod(j+redblack,2))
        do i = int(1)+is,int(3),2
          k2 = F1o8*Volume(i,j)*(Psi(i,j)**(-8))*KK(i,j)
          if((i.eq.int(1)).and.(imin.ne.int(1))) then
            numpsi = LapX(i,j)  *(PsiD(i+1,j) - PsiD(i,j)) +
     .               LapY(i,j)  *(PsiD(i,j+1) - PsiD(i,j)) +
     .               LapY(i,j-1)*(PsiD(i,j-1) - PsiD(i,j)) +
     .               Psi(i,j)*k2
            denpsi = ONE/(LapX(i,j) +
     .                    LapY(i,j) + LapY(i,j-1) +
     .                    SEVEN*k2)
          else if((i.eq.int(3)).and.(imax.ne.int(3))) then
            numpsi = LapX(i-1,j)*(PsiD(i-1,j) - PsiD(i,j)) +
     .               LapY(i,j)  *(PsiD(i,j+1) - PsiD(i,j)) +
     .               LapY(i,j-1)*(PsiD(i,j-1) - PsiD(i,j)) +
     .               Psi(i,j)*k2
            denpsi = ONE/(            LapX(i-1,j) +
     .                    LapY(i,j) + LapY(i,j-1) +
     .                    SEVEN*k2)
          else
            numpsi = LapX(i,j)  *(PsiD(i+1,j) - PsiD(i,j)) +
     .               LapX(i-1,j)*(PsiD(i-1,j) - PsiD(i,j)) +
     .               LapY(i,j)  *(PsiD(i,j+1) - PsiD(i,j)) +
     .               LapY(i,j-1)*(PsiD(i,j-1) - PsiD(i,j)) +
     .               Psi(i,j)*k2
            denpsi = ONE/(LapX(i,j) + LapX(i-1,j) +
     .                    LapY(i,j) + LapY(i,j-1) +
     .                    SEVEN*k2)
          end if
          PsiD(i,j) = PsiD(i,j) + (numpsi + PsiTau(i,j))*denpsi
          Psi(i,j) = PsiD(i,j) + ONE
        end do
      end do
c- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c- Outer Boundary
      if(jmax.ne.int(4)) then
        j = jmax+1
        is = ABS(mod(j+redblack,2))
        do i = int(1)+is,int(3),2
          k2 = F1o8*Volume(i,j)*(Psi(i,j)**(-8))*KK(i,j)
          if((i.eq.int(1)).and.(imin.ne.int(1))) then
            numpsi = LapX(i,j)  *(PsiD(i+1,j) - PsiD(i,j)) -
     .               LapY(i,j)  *(              PsiD(i,j)) +
     .               LapY(i,j-1)*(PsiD(i,j-1) - PsiD(i,j)) +
     .               Psi(i,j)*k2
            denpsi = ONE/(LapX(i,j) +
     .                    LapY(i,j) + LapY(i,j-1) +
     .                    SEVEN*k2)
          else if((i.eq.int(3)).and.(imax.ne.int(3))) then
            numpsi = LapX(i-1,j)*(PsiD(i-1,j) - PsiD(i,j)) -
     .               LapY(i,j)  *(              PsiD(i,j)) +
     .               LapY(i,j-1)*(PsiD(i,j-1) - PsiD(i,j)) +
     .               Psi(i,j)*k2
            denpsi = ONE/(            LapX(i-1,j) +
     .                    LapY(i,j) + LapY(i,j-1) +
     .                    SEVEN*k2)
          else
            numpsi = LapX(i,j)  *(PsiD(i+1,j) - PsiD(i,j)) +
     .               LapX(i-1,j)*(PsiD(i-1,j) - PsiD(i,j)) -
     .               LapY(i,j)  *(              PsiD(i,j)) +
     .               LapY(i,j-1)*(PsiD(i,j-1) - PsiD(i,j)) +
     .               Psi(i,j)*k2
            denpsi = ONE/(LapX(i,j) + LapX(i-1,j) +
     .                    LapY(i,j) + LapY(i,j-1) +
     .                    SEVEN*k2      )
          end if
          PsiD(i,j) = PsiD(i,j) + (numpsi + PsiTau(i,j))*denpsi
          Psi(i,j) = PsiD(i,j) + ONE
        end do
      end if
c
      return
      end
