c#######################################################################
c#
c#                Cadez Smoothers
c#
c#######################################################################
c
c#######################################################################
c# Psi Red/Black Smoother
c#######################################################################
      subroutine psismoothrb(Psi,PsiD,KK,
     .                       PsiTau,LapW,LapX,LapY,Volume,
     .                       ext,int,fil,redblack)
      implicit none
c-----------------------------------------------------------------------
      integer ext(6),int(6),fil(6),redblack
      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 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,is
      real*8 numpsi, denpsi, k2, tmp1,tmp2,bet,ac
      integer NPHIMAX
      parameter(NPHIMAX = 512)
      real*8 ta(2,NPHIMAX)
      integer red,black
      parameter(red = 0,black = 1)
      real*8 F1o8,ONE,TWO,SEVEN
      parameter(F1o8 = 1.D0/8.D0,ONE = 1.D0,TWO = 2.D0,SEVEN = 7.D0)
c
c-----------------------------------------------------------------------
      if(int(4) - int(1) + 1 .gt.NPHIMAX) then
        write(*,*) "psismoothrb:  NPHIMAX too small!"
        stop
      end if
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- Perform Red/Black checkerboard smoothing for Psi
c- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c- Inner Boundary
      if(jmin.ne.int(3)) then
        j = jmin-1
        is = ABS(mod(j+redblack,2))
        do i = int(2)+is,int(5),2
          tmp1 = F1o8*Volume(i,j)
          ac = -LapW(i,j)
          if((i.eq.int(2)).and.(imin.ne.int(2))) then
            tmp2 = TWO*LapW(i,j) +
     .                 LapX(i,j) +
     .                 LapY(i,j) + LapX(i,j-1)*LapY(i,j-1)
            do k = kmin,kmax
              k2 = tmp1*(Psi(k,i,j)**(-8))*KK(k,i,j)
              numpsi = LapW(i,j)  *(PsiD(k+1,i,j) +
     .                              PsiD(k-1,i,j) - TWO*PsiD(k,i,j)) +
     .                 LapX(i,j)  *(PsiD(k,i+1,j) - PsiD(k,i,j)) +
     .                 LapY(i,j)  *(PsiD(k,i,j+1) - PsiD(k,i,j)) -
     .                 LapX(i,j-1)*
     .                 LapY(i,j-1)*(ONE           + PsiD(k,i,j)) +
     .                 Psi(k,i,j)*k2
              ta(2,k-kmin+1) = numpsi + PsiTau(k,i,j)
              ta(1,k-kmin+1) = tmp2 + SEVEN*k2
            end do
            bet = ONE/ta(1,1)
            ta(2,1) = ta(2,1)*bet
            do k = 2,kmax-kmin+1
              ta(1,k-1) = ac*bet
              bet = ONE/(ta(1,k) - ac*ta(1,k-1))
              ta(2,k) = (ta(2,k) - ac*ta(2,k-1))*bet
            end do
            do k = kmax-kmin,1,-1
              ta(2,k) = ta(2,k) - ta(2,k+1)*ta(1,k)
            end do
            do k = kmin,kmax
              PsiD(k,i,j) = PsiD(k,i,j) + ta(2,k-kmin+1)
              Psi(k,i,j) = PsiD(k,i,j) + ONE
            end do
          else if((i.eq.int(5)).and.(imax.ne.int(5))) then
            tmp2 = TWO*LapW(i,j) +
     .                               LapX(i-1,j) +
     .                 LapY(i,j)   + LapX(i,j-1)*LapY(i,j-1)
            do k = kmin,kmax
              k2 = tmp1*(Psi(k,i,j)**(-8))*KK(k,i,j)
              numpsi = LapW(i,j)  *(PsiD(k+1,i,j) +
     .                              PsiD(k-1,i,j) - TWO*PsiD(k,i,j)) +
     .                 LapX(i-1,j)*(PsiD(k,i-1,j) - PsiD(k,i,j)) +
     .                 LapY(i,j)  *(PsiD(k,i,j+1) - PsiD(k,i,j)) -
     .                 LapX(i,j-1)*
     .                 LapY(i,j-1)*(ONE           + PsiD(k,i,j)) +
     .                 Psi(k,i,j)*k2
              ta(2,k-kmin+1) = numpsi + PsiTau(k,i,j)
              ta(1,k-kmin+1) = tmp2 + SEVEN*k2
            end do
            bet = ONE/ta(1,1)
            ta(2,1) = ta(2,1)*bet
            do k = 2,kmax-kmin+1
              ta(1,k-1) = ac*bet
              bet = ONE/(ta(1,k) - ac*ta(1,k-1))
              ta(2,k) = (ta(2,k) - ac*ta(2,k-1))*bet
            end do
            do k = kmax-kmin,1,-1
              ta(2,k) = ta(2,k) - ta(2,k+1)*ta(1,k)
            end do
            do k = kmin,kmax
              PsiD(k,i,j) = PsiD(k,i,j) + ta(2,k-kmin+1)
              Psi(k,i,j) = PsiD(k,i,j) + ONE
            end do
          else
            tmp2 = TWO*LapW(i,j) +
     .                 LapX(i,j) + LapX(i-1,j) +
     .                 LapY(i,j) + LapX(i,j-1)*LapY(i,j-1)
            do k = kmin,kmax
              k2 = tmp1*(Psi(k,i,j)**(-8))*KK(k,i,j)
              numpsi = LapW(i,j)  *(PsiD(k+1,i,j) +
     .                              PsiD(k-1,i,j) - TWO*PsiD(k,i,j)) +
     .                 LapX(i,j)  *(PsiD(k,i+1,j) - PsiD(k,i,j)) +
     .                 LapX(i-1,j)*(PsiD(k,i-1,j) - PsiD(k,i,j)) +
     .                 LapY(i,j)  *(PsiD(k,i,j+1) - PsiD(k,i,j)) -
     .                 LapX(i,j-1)*
     .                 LapY(i,j-1)*(ONE           + PsiD(k,i,j)) +
     .                 Psi(k,i,j)*k2
              ta(2,k-kmin+1) = numpsi + PsiTau(k,i,j)
              ta(1,k-kmin+1) = tmp2 + SEVEN*k2
            end do
            bet = ONE/ta(1,1)
            ta(2,1) = ta(2,1)*bet
            do k = 2,kmax-kmin+1
              ta(1,k-1) = ac*bet
              bet = ONE/(ta(1,k) - ac*ta(1,k-1))
              ta(2,k) = (ta(2,k) - ac*ta(2,k-1))*bet
            end do
            do k = kmax-kmin,1,-1
              ta(2,k) = ta(2,k) - ta(2,k+1)*ta(1,k)
            end do
            do k = kmin,kmax
              PsiD(k,i,j) = PsiD(k,i,j) + ta(2,k-kmin+1)
              Psi(k,i,j) = PsiD(k,i,j) + ONE
            end do
          end if
        end do
      end if
c- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c- Interior
      do j = jmin,jmax
        is = ABS(mod(j+redblack,2))
        do i = int(2)+is,int(5),2
          tmp1 = F1o8*Volume(i,j)
          ac = -LapW(i,j)
          if((i.eq.int(2)).and.(imin.ne.int(2))) then
            tmp2 = TWO*LapW(i,j) +
     .                 LapX(i,j) +
     .                 LapY(i,j) + LapY(i,j-1)
            do k = kmin,kmax
              k2 = tmp1*(Psi(k,i,j)**(-8))*KK(k,i,j)
              numpsi = LapW(i,j)  *(PsiD(k+1,i,j) +
     .                              PsiD(k-1,i,j) - TWO*PsiD(k,i,j)) +
     .                 LapX(i,j)  *(PsiD(k,i+1,j) - PsiD(k,i,j)) +
     .                 LapY(i,j)  *(PsiD(k,i,j+1) - PsiD(k,i,j)) +
     .                 LapY(i,j-1)*(PsiD(k,i,j-1) - PsiD(k,i,j)) +
     .                 Psi(k,i,j)*k2
              ta(2,k-kmin+1) = numpsi + PsiTau(k,i,j)
              ta(1,k-kmin+1) = tmp2 + SEVEN*k2
            end do
            bet = ONE/ta(1,1)
            ta(2,1) = ta(2,1)*bet
            do k = 2,kmax-kmin+1
              ta(1,k-1) = ac*bet
              bet = ONE/(ta(1,k) - ac*ta(1,k-1))
              ta(2,k) = (ta(2,k) - ac*ta(2,k-1))*bet
            end do
            do k = kmax-kmin,1,-1
              ta(2,k) = ta(2,k) - ta(2,k+1)*ta(1,k)
            end do
            do k = kmin,kmax
              PsiD(k,i,j) = PsiD(k,i,j) + ta(2,k-kmin+1)
              Psi(k,i,j) = PsiD(k,i,j) + ONE
            end do
          else if((i.eq.int(5)).and.(imax.ne.int(5))) then
            tmp2 = TWO*LapW(i,j) +
     .                             LapX(i-1,j) +
     .                 LapY(i,j) + LapY(i,j-1)
            do k = kmin,kmax
              k2 = tmp1*(Psi(k,i,j)**(-8))*KK(k,i,j)
              numpsi = LapW(i,j)  *(PsiD(k+1,i,j) +
     .                              PsiD(k-1,i,j) - TWO*PsiD(k,i,j)) +
     .                 LapX(i-1,j)*(PsiD(k,i-1,j) - PsiD(k,i,j)) +
     .                 LapY(i,j)  *(PsiD(k,i,j+1) - PsiD(k,i,j)) +
     .                 LapY(i,j-1)*(PsiD(k,i,j-1) - PsiD(k,i,j)) +
     .                 Psi(k,i,j)*k2
              ta(2,k-kmin+1) = numpsi + PsiTau(k,i,j)
              ta(1,k-kmin+1) = tmp2 + SEVEN*k2
            end do
            bet = ONE/ta(1,1)
            ta(2,1) = ta(2,1)*bet
            do k = 2,kmax-kmin+1
              ta(1,k-1) = ac*bet
              bet = ONE/(ta(1,k) - ac*ta(1,k-1))
              ta(2,k) = (ta(2,k) - ac*ta(2,k-1))*bet
            end do
            do k = kmax-kmin,1,-1
              ta(2,k) = ta(2,k) - ta(2,k+1)*ta(1,k)
            end do
            do k = kmin,kmax
              PsiD(k,i,j) = PsiD(k,i,j) + ta(2,k-kmin+1)
              Psi(k,i,j) = PsiD(k,i,j) + ONE
            end do
          else
            tmp2 = TWO*LapW(i,j) +
     .                 LapX(i,j) + LapX(i-1,j) +
     .                 LapY(i,j) + LapY(i,j-1)
            do k = kmin,kmax
              k2 = tmp1*(Psi(k,i,j)**(-8))*KK(k,i,j)
              numpsi = LapW(i,j)  *(PsiD(k+1,i,j) +
     .                              PsiD(k-1,i,j) - TWO*PsiD(k,i,j)) +
     .                 LapX(i,j)  *(PsiD(k,i+1,j) - PsiD(k,i,j)) +
     .                 LapX(i-1,j)*(PsiD(k,i-1,j) - PsiD(k,i,j)) +
     .                 LapY(i,j)  *(PsiD(k,i,j+1) - PsiD(k,i,j)) +
     .                 LapY(i,j-1)*(PsiD(k,i,j-1) - PsiD(k,i,j)) +
     .                 Psi(k,i,j)*k2
              ta(2,k-kmin+1) = numpsi + PsiTau(k,i,j)
              ta(1,k-kmin+1) = tmp2 + SEVEN*k2
            end do
            bet = ONE/ta(1,1)
            ta(2,1) = ta(2,1)*bet
            do k = 2,kmax-kmin+1
              ta(1,k-1) = ac*bet
              bet = ONE/(ta(1,k) - ac*ta(1,k-1))
              ta(2,k) = (ta(2,k) - ac*ta(2,k-1))*bet
            end do
            do k = kmax-kmin,1,-1
              ta(2,k) = ta(2,k) - ta(2,k+1)*ta(1,k)
            end do
            do k = kmin,kmax
              PsiD(k,i,j) = PsiD(k,i,j) + ta(2,k-kmin+1)
              Psi(k,i,j) = PsiD(k,i,j) + ONE
            end do
          end if
        end do
      end do
c- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c- Outer Boundary
      if(jmax.ne.int(6)) then
        j = jmax+1
        is = ABS(mod(j+redblack,2))
        do i = int(2)+is,int(5),2
          tmp1 = F1o8*Volume(i,j)
          ac = -LapW(i,j)
          if((i.eq.int(2)).and.(imin.ne.int(2))) then
            tmp2 = TWO*LapW(i,j) +
     .                 LapX(i,j) +
     .                 LapY(i,j) + LapY(i,j-1)
            do k = kmin,kmax
              k2 = tmp1*(Psi(k,i,j)**(-8))*KK(k,i,j)
              numpsi = LapW(i,j)  *(PsiD(k+1,i,j) +
     .                              PsiD(k-1,i,j) - TWO*PsiD(k,i,j)) +
     .                 LapX(i,j)  *(PsiD(k,i+1,j) - PsiD(k,i,j)) -
     .                 LapY(i,j)  *(                PsiD(k,i,j)) +
     .                 LapY(i,j-1)*(PsiD(k,i,j-1) - PsiD(k,i,j)) +
     .                 Psi(k,i,j)*k2
              ta(2,k-kmin+1) = numpsi + PsiTau(k,i,j)
              ta(1,k-kmin+1) = tmp2 + SEVEN*k2
            end do
            bet = ONE/ta(1,1)
            ta(2,1) = ta(2,1)*bet
            do k = 2,kmax-kmin+1
              ta(1,k-1) = ac*bet
              bet = ONE/(ta(1,k) - ac*ta(1,k-1))
              ta(2,k) = (ta(2,k) - ac*ta(2,k-1))*bet
            end do
            do k = kmax-kmin,1,-1
              ta(2,k) = ta(2,k) - ta(2,k+1)*ta(1,k)
            end do
            do k = kmin,kmax
              PsiD(k,i,j) = PsiD(k,i,j) + ta(2,k-kmin+1)
              Psi(k,i,j) = PsiD(k,i,j) + ONE
            end do
          else if((i.eq.int(5)).and.(imax.ne.int(5))) then
            tmp2 = TWO*LapW(i,j) +
     .                             LapX(i-1,j) +
     .                 LapY(i,j) + LapY(i,j-1)
            do k = kmin,kmax
              k2 = tmp1*(Psi(k,i,j)**(-8))*KK(k,i,j)
              numpsi = LapW(i,j)  *(PsiD(k+1,i,j) +
     .                              PsiD(k-1,i,j) - TWO*PsiD(k,i,j)) +
     .                 LapX(i-1,j)*(PsiD(k,i-1,j) - PsiD(k,i,j)) -
     .                 LapY(i,j)  *(                PsiD(k,i,j)) +
     .                 LapY(i,j-1)*(PsiD(k,i,j-1) - PsiD(k,i,j)) +
     .                 Psi(k,i,j)*k2
              ta(2,k-kmin+1) = numpsi + PsiTau(k,i,j)
              ta(1,k-kmin+1) = tmp2 + SEVEN*k2
            end do
            bet = ONE/ta(1,1)
            ta(2,1) = ta(2,1)*bet
            do k = 2,kmax-kmin+1
              ta(1,k-1) = ac*bet
              bet = ONE/(ta(1,k) - ac*ta(1,k-1))
              ta(2,k) = (ta(2,k) - ac*ta(2,k-1))*bet
            end do
            do k = kmax-kmin,1,-1
              ta(2,k) = ta(2,k) - ta(2,k+1)*ta(1,k)
            end do
            do k = kmin,kmax
              PsiD(k,i,j) = PsiD(k,i,j) + ta(2,k-kmin+1)
              Psi(k,i,j) = PsiD(k,i,j) + ONE
            end do
          else
            tmp2 = TWO*LapW(i,j) +
     .                 LapX(i,j) + LapX(i-1,j) +
     .                 LapY(i,j) + LapY(i,j-1)
            do k = kmin,kmax
              k2 = tmp1*(Psi(k,i,j)**(-8))*KK(k,i,j)
              numpsi = LapW(i,j)  *(PsiD(k+1,i,j) +
     .                              PsiD(k-1,i,j) - TWO*PsiD(k,i,j)) +
     .                 LapX(i,j)  *(PsiD(k,i+1,j) - PsiD(k,i,j)) +
     .                 LapX(i-1,j)*(PsiD(k,i-1,j) - PsiD(k,i,j)) -
     .                 LapY(i,j)  *(                PsiD(k,i,j)) +
     .                 LapY(i,j-1)*(PsiD(k,i,j-1) - PsiD(k,i,j)) +
     .                 Psi(k,i,j)*k2
              ta(2,k-kmin+1) = numpsi + PsiTau(k,i,j)
              ta(1,k-kmin+1) = tmp2 + SEVEN*k2
            end do
            bet = ONE/ta(1,1)
            ta(2,1) = ta(2,1)*bet
            do k = 2,kmax-kmin+1
              ta(1,k-1) = ac*bet
              bet = ONE/(ta(1,k) - ac*ta(1,k-1))
              ta(2,k) = (ta(2,k) - ac*ta(2,k-1))*bet
            end do
            do k = kmax-kmin,1,-1
              ta(2,k) = ta(2,k) - ta(2,k+1)*ta(1,k)
            end do
            do k = kmin,kmax
              PsiD(k,i,j) = PsiD(k,i,j) + ta(2,k-kmin+1)
              Psi(k,i,j) = PsiD(k,i,j) + ONE
            end do
          end if
        end do
      end if
c
      return
      end
