      subroutine psi_ds(PsiD,Psi,KK,
     .                  PsiRes,LapX,LapY,Volume,A,B,ext,int,KD)
      implicit none
      integer ext(4),int(4),KD
      real*8 PsiD(ext(1):ext(3),ext(2):ext(4))
      real*8 Psi(ext(1):ext(3),ext(2):ext(4))
      real*8 KK(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))
      real*8 A(KD+1,int(1):int(3),int(2):int(4))
      real*8 B(int(1):int(3),int(2):int(4))
      integer i,j,k,imin,imax,jmin,jmax,N,info,jmina,jmaxa
      real*8 k2
      real*8 ZERO,F7o8,THREE
      parameter (ZERO = 0.D0,F7o8 = 7.D0/8.D0,THREE = 3.D0)
c
      imin = int(1)
      jmin = int(2)
      imax = int(3)
      jmax = int(4)
      N = (jmax - jmin + 1)*(imax - imin + 1)
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
        k2 = F7o8*(Psi(i,j)**(-8))*Volume(i,j)*KK(i,j)
        A(1,i,j) = LapX(i,j) + LapX(i,j-1)*LapY(i,j-1) + LapY(i,j)
     .           + k2
        A(2,i,j) = -LapX(i,j)
        do k=3,KD
          A(k,i,j) = ZERO
        end do
        A(KD+1,i,j) = -LapY(i,j)
        do i=imin+1,imax-1
          k2 = F7o8*(Psi(i,j)**(-8))*Volume(i,j)*KK(i,j)
          A(1,i,j) = LapX(i-1,j) + LapX(i,j)
     .             + LapX(i,j-1)*LapY(i,j-1) + LapY(i,j)
     .             + k2
          A(2,i,j) = -LapX(i,j)
          do k=3,KD
            A(k,i,j) = ZERO
          end do
          A(KD+1,i,j) = -LapY(i,j)
        end do
        i = imax
        k2 = F7o8*(Psi(i,j)**(-8))*Volume(i,j)*KK(i,j)
        A(1,i,j) = LapX(i-1,j) + LapX(i,j-1)*LapY(i,j-1) + LapY(i,j)
     .           + k2
        do k=2,KD
          A(k,i,j) = ZERO
        end do
        A(KD+1,i,j) = -LapY(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
        k2 = F7o8*(Psi(i,j)**(-8))*Volume(i,j)*KK(i,j)
        A(1,i,j) = LapX(i,j) + LapY(i,j-1) + LapY(i,j)
     .           + k2
        A(2,i,j) = -LapX(i,j)
        do k=3,KD+1
          A(k,i,j) = ZERO
        end do
        do i=imin+1,imax-1
          k2 = F7o8*(Psi(i,j)**(-8))*Volume(i,j)*KK(i,j)
          A(1,i,j) = LapX(i-1,j) + LapX(i,j)
     .             + LapY(i,j-1) + LapY(i,j)
     .             + k2
          A(2,i,j) = -LapX(i,j)
          do k=3,KD+1
            A(k,i,j) = ZERO
          end do
        end do
        i = imax
        k2 = F7o8*(Psi(i,j)**(-8))*Volume(i,j)*KK(i,j)
        A(1,i,j) = LapX(i-1,j) + LapY(i,j-1) + LapY(i,j)
     .           + k2
        do k=2,KD+1
          A(k,i,j) = ZERO
        end do
        jmaxa = jmax-1
      endif
c- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c- Interior
      do j=jmina,jmaxa
        i = imin
        k2 = F7o8*(Psi(i,j)**(-8))*Volume(i,j)*KK(i,j)
        A(1,i,j) = LapX(i,j) + LapY(i,j-1) + LapY(i,j)
     .           + k2
        A(2,i,j) = -LapX(i,j)
        do k=3,KD
          A(k,i,j) = ZERO
        end do
        A(KD+1,i,j) = -LapY(i,j)
        do i=imin+1,imax-1
          k2 = F7o8*(Psi(i,j)**(-8))*Volume(i,j)*KK(i,j)
          A(1,i,j) = LapX(i-1,j) + LapX(i,j)
     .             + LapY(i,j-1) + LapY(i,j)
     .             + k2
          A(2,i,j) = -LapX(i,j)
          do k=3,KD
            A(k,i,j) = ZERO
          end do
          A(KD+1,i,j) = -LapY(i,j)
        end do
        i = imax
        k2 = F7o8*(Psi(i,j)**(-8))*Volume(i,j)*KK(i,j)
        A(1,i,j) = LapX(i-1,j) + LapY(i,j-1) + LapY(i,j)
     .           + k2
        do k=2,KD
          A(k,i,j) = ZERO
        end do
        A(KD+1,i,j) = -LapY(i,j)
      end do
c
c-----------------------------------------------------------------------
c- Fill the RHS
c-----------------------------------------------------------------------
      do j=jmin,jmax
        do i=imin,imax
          B(i,j) = PsiRes(i,j)
        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
          PsiD(i,j) = PsiD(i,j) - B(i,j)
        end do
      end do
c
      return
      end
