!=======================================================================
      subroutine sip5(ilumode,nx,ny,itmax,iterST,level,ierr,
     &             alpha,tol,A,AI,X,B,wksp)
!=======================================================================
      implicit none
      integer ilumode,nx,ny,itmax,iterST,level,ierr
      real*8  alpha,tol
      real*8  A(5,nx,ny),AI(5,nx,ny)
      real*8  X(nx,ny),B(nx,ny),wksp(nx,ny)

!-----------------------------------------------------------------------
! Strongly Implicit Procedure:
! Stone, H.L. (1968), SIAM J. Numer. Anal. 5, pp. 530-558.
!
! A:    the 5-point matrix from structured grid in 2D
! AI:   used for the Stone's ILU matrix
! X:    the initial guess (input) and the unknown (output)
! b:    the RHS of the algebraic equation
! wksp: used as the "residual" and the "correction term".
!-----------------------------------------------------------------------

!---- Local variables

      integer ix,iy,level0
      real*8  l2norm,l8norm,res,res0

      level0=0

!---- Informing: stdout

      if(level.ge.1)then
         print'("SIP5:",$)'
         print*,"nx=",nx," ny=",ny
      endif

!---- Make the ILU matrix, using Stone's method

      if(ilumode.eq.1) then
         call sip_ilu(nx,ny,level0,ierr,alpha,A,AI)
      endif


!-----------------------------
      do iterST=1,itmax
!-----------------------------

!---- Get Residual

      call residual5(nx,ny,A,X,B,wksp)


!---- Compute Correction

      call sip_subst(nx,ny,AI,wksp)


!---- Add the correction to the solution

      do iy=1,ny
      do ix=1,nx
         X(ix,iy)=X(ix,iy)+wksp(ix,iy)
      enddo
      enddo


!---- Check to stop

      res=l8norm(nx*ny,wksp)
      if(iterST.eq.1)then
         res0=res
         if(res0.lt.1.d-10) goto 9999
      endif
      if((res/res0).le.tol) goto 9999

!-----------------------------
      enddo
!-----------------------------

 9999 continue

      return
      end


!=======================================================================
      subroutine sip_ilu(nx,ny,level,ierr,alpha,A,AI)
!=======================================================================
      implicit none
      integer nx,ny,level,ierr
      real*8  alpha
      real*8  A(5,nx,ny),AI(5,nx,ny)

!---- Local variables

      integer ix,iy
      real*8  LS,LW,LP,UE,UN
      real*8  UE1,UE2,UN1,UN2  ! UE1=UE(ix-1,iy); UE2=UE(ix,iy-1)
      real*8  zero,one

!---- Local setting

      zero=0.d0
      one =1.d0

!---- Informing: stdout

      if(level.ge.1)then
         print'("SIP_ILU:",$)'
         print*,"nx=",nx," ny=",ny
      endif

!---- Make the ILU matrix, using Stone's method

      do iy=1,ny
      do ix=1,nx

         if(ix.eq.1)then
            UE1=zero
         else
            UE1=AI(4,ix-1,iy)
         endif

         if(iy.eq.1)then
            UE2=zero
         else
            UE2=AI(4,ix,iy-1)
         endif

         if(ix.eq.1)then
            UN1=zero
         else
            UN1=AI(5,ix-1,iy)
         endif

         if(iy.eq.1)then
            UN2=zero
         else
            UN2=AI(5,ix,iy-1)
         endif

         LS=A(1,ix,iy)/(one+alpha*UE2)
         LW=A(2,ix,iy)/(one+alpha*UN1)
         LP=A(3,ix,iy)+alpha*(LS*UE2+LW*UN1)-(LS*UN2+LW*UE1)
         UE=(A(4,ix,iy)-alpha*LS*UE2)/LP
         UN=(A(5,ix,iy)-alpha*LW*UN1)/LP

         AI(1,ix,iy)=LS
         AI(2,ix,iy)=LW
         AI(3,ix,iy)=LP
         AI(4,ix,iy)=UE
         AI(5,ix,iy)=UN

      enddo
      enddo

      return
      end


!=======================================================================
      subroutine sip_subst(nx,ny,ailu,x)
!=======================================================================
      implicit none
      integer nx,ny
      real*8  ailu(5,nx,ny),x(nx,ny)

!---- Here the main diagonal of "U" is one.

!---- Local variables

      integer ix,iy

!---- Forward substitution

      iy=1
         do ix=2,nx
            x(ix,iy)=(x(ix,iy)-ailu(2,ix,iy)*x(ix-1,iy))/ailu(3,ix,iy)
         enddo

      do iy=2,ny
         ix=1
            x(ix,iy)=(x(ix,iy)-ailu(1,ix,iy)*x(ix,iy-1))/ailu(3,ix,iy)
         do ix=2,nx
            x(ix,iy)=(x(ix,iy)-ailu(1,ix,iy)*x(ix,iy-1)
     &                        -ailu(2,ix,iy)*x(ix-1,iy))/ailu(3,ix,iy)
         enddo
      enddo

!---- Backward substitution

      iy=ny
         do ix=nx-1,1,-1
            x(ix,iy)=x(ix,iy)-ailu(4,ix,iy)*x(ix+1,iy)
         enddo

      do iy=ny-1,1,-1
         ix=nx
            x(ix,iy)=x(ix,iy)-ailu(5,ix,iy)*x(ix,iy+1)
         do ix=nx-1,1,-1
            x(ix,iy)=x(ix,iy)-ailu(4,ix,iy)*x(ix+1,iy)
     &                       -ailu(5,ix,iy)*x(ix,iy+1)
         enddo
      enddo

      return
      end


c ======================================================================
      subroutine residual5(nx,ny,a,x,f,res)
c ======================================================================
      implicit none
      integer nx,ny
      real*8  a(5,nx,ny),x(nx,ny),f(nx,ny),res(nx,ny)

      integer i,j,jm1,jp1
      integer nxm1,nym1
      real*8  axij

c---- res = f - a x

      nxm1=nx-1
      nym1=ny-1

c---- interior points

      do j=2,nym1
         jm1=j-1
         jp1=j+1
      do i=2,nxm1
         axij=a(1,i,j)*x(i,jm1)+a(2,i,j)*x(i-1,j)+a(3,i,j)*x(i,j)
     &       +a(4,i,j)*x(i+1,j)+a(5,i,j)*x(i,jp1)
         res(i,j)=f(i,j)-axij
      end do
      end do

c---- four corners

      i=1
      j=1
      axij=a(3,i,j)*x(i,j)+a(4,i,j)*x(i+1,j)+a(5,i,j)*x(i,j+1)
      res(i,j)=f(i,j)-axij

      i=nx
      j=1
      axij=a(3,i,j)*x(i,j)+a(2,i,j)*x(i-1,j)+a(5,i,j)*x(i,j+1)
      res(i,j)=f(i,j)-axij

      i=1
      j=ny
      axij=a(3,i,j)*x(i,j)+a(1,i,j)*x(i,j-1)+a(4,i,j)*x(i+1,j)
      res(i,j)=f(i,j)-axij

      i=nx
      j=ny
      axij=a(3,i,j)*x(i,j)+a(1,i,j)*x(i,j-1)+a(2,i,j)*x(i-1,j)
      res(i,j)=f(i,j)-axij

c---- four sides

      i=1
      do j=2,nym1
         axij=a(1,i,j)*x(i,j-1)+a(3,i,j)*x(i,j)
     &       +a(4,i,j)*x(i+1,j)+a(5,i,j)*x(i,j+1)
         res(i,j)=f(i,j)-axij
      end do

      i=nx
      do j=2,nym1
         axij=a(1,i,j)*x(i,j-1)+a(2,i,j)*x(i-1,j)+a(3,i,j)*x(i,j)
     &       +a(5,i,j)*x(i,j+1)
         res(i,j)=f(i,j)-axij
      end do

      j=1
      do i=2,nxm1
         axij=a(2,i,j)*x(i-1,j)+a(3,i,j)*x(i,j)
     &       +a(4,i,j)*x(i+1,j)+a(5,i,j)*x(i,j+1)
         res(i,j)=f(i,j)-axij
      end do

      j=ny
      do i=2,nxm1
         axij=a(1,i,j)*x(i,j-1)+a(2,i,j)*x(i-1,j)+a(3,i,j)*x(i,j)
     &       +a(4,i,j)*x(i+1,j)
         res(i,j)=f(i,j)-axij
      end do


      return
      end