      subroutine psi_ds(PsiD,Psi,KK,
     .                  PsiRes,LapW,LapX,LapY,Volume,A,B,ext,int,KD)
      implicit none
      integer ext(6),int(6),KD
      real*8 PsiD(ext(1):ext(4),ext(2):ext(5),ext(3):ext(6))
      real*8 Psi(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 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))
      real*8 A(KD+1,int(1):int(4),int(2):int(5),int(3):int(6))
      real*8 B(int(1):int(4),int(2):int(5),int(3):int(6))
      integer k,i,j,l,kmin,kmax,imin,imax,jmin,jmax
      integer Nk,N,info,jmina,jmaxa
      real*8 k2, tmp1, tmp2
      real*8 ZERO,F7o8,TWO,THREE
      parameter (ZERO = 0.D0,F7o8 = 7.D0/8.D0,TWO = 2.D0,THREE = 3.D0)
c
      kmin = int(1)
      imin = int(2)
      jmin = int(3)
      kmax = int(4)
      imax = int(5)
      jmax = int(6)
      Nk = kmax - kmin + 1
      N = (jmax - jmin + 1)*(imax - imin + 1)*Nk
c
c-----------------------------------------------------------------------
c- Fill A
c-----------------------------------------------------------------------
c- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c- Inner boundary !! Handle inner boundary as special !!
c- (in Cadez inner boundary has jmin < 1)
      jmina = jmin
      if(jmin.lt.1) then
        j = jmin
        i = imin
        tmp1 = TWO*LapW(i,j)
     .       + LapX(i,j) + LapX(i,j-1)*LapY(i,j-1) + LapY(i,j)
        tmp2 = F7o8*Volume(i,j)
        do k=kmin,kmax
          A(1,k,i,j) = tmp1 + tmp2*(Psi(k,i,j)**(-8))*KK(k,i,j)
          A(2,k,i,j) = -LapW(i,j)
          do l=3,Nk
            A(l,k,i,j) = ZERO
          end do
          A(Nk+1,k,i,j) = -LapX(i,j)
          do l=Nk+2,KD
            A(l,k,i,j) = ZERO
          end do
          A(KD+1,k,i,j) = -LapY(i,j)
        end do
        A(2,kmax,i,j) = ZERO
        A(Nk,kmin,i,j) = -LapW(i,j)
        do i=imin+1,imax-1
          tmp1 = TWO*LapW(i,j)
     .         + LapX(i-1,j) + LapX(i,j)
     .         + LapX(i,j-1)*LapY(i,j-1) + LapY(i,j)
          tmp2 = F7o8*Volume(i,j)
          do k=kmin,kmax
            A(1,k,i,j) = tmp1 + tmp2*(Psi(k,i,j)**(-8))*KK(k,i,j)
            A(2,k,i,j) = -LapW(i,j)
            do l=3,Nk
              A(l,k,i,j) = ZERO
            end do
            A(Nk+1,k,i,j) = -LapX(i,j)
            do l=Nk+2,KD
              A(l,k,i,j) = ZERO
            end do
            A(KD+1,k,i,j) = -LapY(i,j)
          end do
          A(2,kmax,i,j) = ZERO
          A(Nk,kmin,i,j) = -LapW(i,j)
        end do
        i = imax
        tmp1 = TWO*LapW(i,j)
     .       + LapX(i-1,j) + LapX(i,j-1)*LapY(i,j-1) + LapY(i,j)
        tmp2 = F7o8*Volume(i,j)
        do k=kmin,kmax
          A(1,k,i,j) = tmp1 + tmp2*(Psi(k,i,j)**(-8))*KK(k,i,j)
          A(2,k,i,j) = -LapW(i,j)
          do l=3,KD
            A(l,k,i,j) = ZERO
          end do
          A(KD+1,k,i,j) = -LapY(i,j)
        end do
        A(2,kmax,i,j) = ZERO
        A(Nk,kmin,i,j) = -LapW(i,j)
        jmina = jmin+1
      end if
c- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c- Outer boundary !! Handle outer boundary as special !!
c- (in Cadez outer boundary has jmax > 0)
      jmaxa = jmax
      if(jmax.gt.0) then
        j = jmax
        i = imin
        tmp1 = TWO*LapW(i,j)
     .       + LapX(i,j) + LapY(i,j-1) + LapY(i,j)
        tmp2 = F7o8*Volume(i,j)
        do k=kmin,kmax
          A(1,k,i,j) = tmp1 + tmp2*(Psi(k,i,j)**(-8))*KK(k,i,j)
          A(2,k,i,j) = -LapW(i,j)
          do l=3,Nk
            A(l,k,i,j) = ZERO
          end do
          A(Nk+1,k,i,j) = -LapX(i,j)
          do l=Nk+2,KD+1
            A(l,k,i,j) = ZERO
          end do
        end do
        A(2,kmax,i,j) = ZERO
        A(Nk,kmin,i,j) = -LapW(i,j)
        do i=imin+1,imax-1
          tmp1 = TWO*LapW(i,j)
     .         + LapX(i-1,j) + LapX(i,j)
     .         + LapY(i,j-1) + LapY(i,j)
          tmp2 = F7o8*Volume(i,j)
          do k=kmin,kmax
            A(1,k,i,j) = tmp1 + tmp2*(Psi(k,i,j)**(-8))*KK(k,i,j)
            A(2,k,i,j) = -LapW(i,j)
            do l=3,Nk
              A(l,k,i,j) = ZERO
            end do
            A(Nk+1,k,i,j) = -LapX(i,j)
            do l=Nk+2,KD+1
              A(l,k,i,j) = ZERO
            end do
          end do
          A(2,kmax,i,j) = ZERO
          A(Nk,kmin,i,j) = -LapW(i,j)
        end do
        i = imax
        tmp1 = TWO*LapW(i,j)
     .       + LapX(i-1,j) + LapY(i,j-1) + LapY(i,j)
        tmp2 = F7o8*Volume(i,j)
        do k=kmin,kmax
          A(1,k,i,j) = tmp1 + tmp2*(Psi(k,i,j)**(-8))*KK(k,i,j)
          A(2,k,i,j) = -LapW(i,j)
          do l=3,KD+1
            A(l,k,i,j) = ZERO
          end do
        end do
        A(2,kmax,i,j) = ZERO
        A(Nk,kmin,i,j) = -LapW(i,j)
        jmaxa = jmax-1
      endif
c- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c- Interior
      do j=jmina,jmaxa
        i = imin
        tmp1 = TWO*LapW(i,j)
     .           + LapX(i,j) + LapY(i,j-1) + LapY(i,j)
        tmp2 = F7o8*Volume(i,j)
        do k=kmin,kmax
          A(1,k,i,j) = tmp1 + tmp2*(Psi(k,i,j)**(-8))*KK(k,i,j)
          A(2,k,i,j) = -LapW(i,j)
          do l=3,Nk
            A(l,k,i,j) = ZERO
          end do
          A(Nk+1,k,i,j) = -LapX(i,j)
          do l=Nk+2,KD
            A(l,k,i,j) = ZERO
          end do
          A(KD+1,k,i,j) = -LapY(i,j)
        end do
        A(2,kmax,i,j) = ZERO
        A(Nk,kmin,i,j) = -LapW(i,j)
        do i=imin+1,imax-1
          tmp1 = TWO*LapW(i,j)
     .         + LapX(i-1,j) + LapX(i,j)
     .         + LapY(i,j-1) + LapY(i,j)
          tmp2 = F7o8*Volume(i,j)
          do k=kmin,kmax
            A(1,k,i,j) = tmp1 + tmp2*(Psi(k,i,j)**(-8))*KK(k,i,j)
            A(2,k,i,j) = -LapW(i,j)
            do l=3,Nk
              A(l,k,i,j) = ZERO
            end do
            A(Nk+1,k,i,j) = -LapX(i,j)
            do l=Nk+2,KD
              A(l,k,i,j) = ZERO
            end do
            A(KD+1,k,i,j) = -LapY(i,j)
          end do
          A(2,kmax,i,j) = ZERO
          A(Nk,kmin,i,j) = -LapW(i,j)
        end do
        i = imax
        tmp1 = TWO*LapW(i,j)
     .           + LapX(i-1,j) + LapY(i,j-1) + LapY(i,j)
        tmp2 = F7o8*Volume(i,j)
        do k=kmin,kmax
          A(1,k,i,j) = tmp1 + tmp2*(Psi(k,i,j)**(-8))*KK(k,i,j)
          A(2,k,i,j) = -LapW(i,j)
          do l=3,KD
            A(l,k,i,j) = ZERO
          end do
          A(KD+1,k,i,j) = -LapY(i,j)
        end do
        A(2,kmax,i,j) = ZERO
        A(Nk,kmin,i,j) = -LapW(i,j)
      end do
c
c-----------------------------------------------------------------------
c- Fill the RHS
c-----------------------------------------------------------------------
      do j=jmin,jmax
        do i=imin,imax
          do k=kmin,kmax
            B(k,i,j) = PsiRes(k,i,j)
          end do
        end do
      end do
c
c-----------------------------------------------------------------------
c- Solve for the correction to Psi
c-----------------------------------------------------------------------
      call dpbsv('L',N,KD,1,A,KD+1,B,N,info)
      if(info.ne.0) then
        write(*,*) "Region 1:2:3 DPBSV error, info = ",info
        return
      end if
c
c-----------------------------------------------------------------------
c- Update the value of Psi
c-----------------------------------------------------------------------
      do j=jmin,jmax
        do i=imin,imax
          do k=kmin,kmax
            PsiD(k,i,j) = PsiD(k,i,j) - B(k,i,j)
          end do
        end do
      end do
c
      return
      end
