PROGRAM STREAMF
! Stream function and vorticity solution of Navier-Stokes equation in 2-D
! Specied tangential and normal components of velocity
! s(i,j) is the z component of the vector potentialor the stream function
! w(i,j) is the z component of the vorticity
! u(i,j) and v(i,j) are x and y components of velocity, respectively
! (x,y)=(xi,yj)
! sbx(ib,j) are the values at the x=(0,1) boundaries; (ib)=(1,2)
! sby(i,jb) are the values at the y=(0,1) boundaries; (jb)=(1,2)
! Tangental component of velocity is specified at the boundaries as follows
! vtbx(ib,j) are the values at the x=(0,1) boundaries; (ib)=(1,2)
! vtby(i,jb) are the values at the y=(0,1) boundaries; (jb)=(1,2)
! Normal components of velocity are specified at the boundaries as follows
! vnbx(ib,j) are the values at the x=(0,1) boundaries; (ib)=(1,2)
! vnby(i,jb) are the values at the y=(0,1) boundaries; (jb)=(1,2)
! Note: Integral of the normal component of velocity must be zero
!
! a, b, c, d, e, and f are the coefficients of the difference equations.
!    a*s(i+1,j)+b*s(i-1,j)+c*s(i,j+1)+d*s(i,j-1)+e*s(i,j)=f
!    a*w(i+1,j)+b*w(i-1,j)+c*w(i,j+1)+d*w(i,j-1)+e*w(i,j)=f
! a=a(i,j,m,n), same for b ... f
! (i,j) are grid indices; (m,n) are (equation index, dependent variable index)
! Index 1 is for stream function; index 2 for vorticity
! e.g., a(i,j,1,1) is coefficient in stream function equation for stream function
!       a(i,j,1,2) is coefficient in stream function equation for vorticity
! Equation 1: a(i,j,1,1)*s(i+1,j)+a(i,j,1,2)*w(i+1,j) + ...
! Equation 2: a(i,j,2,1)*s(i+1,j)+a(i,j,2,2)*w(i+1,j) + ...
!
! JMAX is the total number of gridpoints in the X and Y directions, inc. boundaries
! JM2 is the number of interior gridpoints in the X and Y directions.
! (i & j)=(1 & JM1) are the boundaries
! JM1=JMAX-1, JM2=JMAX-2
! del is delta X & Y : del=1/(JMAX-1)
! alf is the aspect ratio, alf=Lx/Ly
! PBCG from Numerical Recipes is used to solve linear system of equations

	INCLUDE 'include2.inc'
	OPEN( 8, FILE='input.dat')
      OPEN(12, FILE='matrix.dat')
      OPEN(13, FILE='u.dat')
      OPEN(14, FILE='v.dat')
      OPEN(19, FILE='s.dat')
      OPEN(20, FILE='w.dat')
      if(NMAX.lt.10*JMAX2)then
        write(*,'(/1x,a)') 'NMAX is too small. Redimension'
        stop
      endif
!
	CALL bc2
	CALL coef2
!     Store elements of matrix as row-indexed sparse storage mode
!     IA is the row index of coefficient matrix, (1,2*JM2*JM2)
!     i & j are index of grid array, (2,JM1)
!     i index is changed most rapidly
!     k is index of row-indexed sparse array
!     ija(k) is the column number of coefficient in coefficient matrix
!     zero elements of near off-diagonal are stored, make sure they are initialized.
!
      NN=2*JM2*JM2
      j=2
      i=1
      do IA=1,NN-1,2
        i=i+1
        if(i.gt.JM1)then
          i=2
          j=j+1
        endif
        do m=1,2
          if(m.eq.1)then
            sa(IA)=e(i,j,m,m)     !  Store diagonal elements
          else
            sa(IA+1)=e(i,j,m,m)   !  Store diagonal elements
          endif
        enddo
      enddo
      ija(1)=NN+2            !  Store size of matrix + 2
      k=NN+1
      j=2
      i=1
      do IA=1,NN-1,2        !  IA is sequence number for first equation of pair
                            !  i.e., It is row counter for coefficient matrix
        i=i+1
        if(i.gt.JM1)then
          i=2
          j=j+1
        endif
!     Store off-diagonal coefficients
        do m=1,2
          k=k+1
          if(m.eq.1)then
            sa(k)=e(i,j,1,2)  ! Store coupling on diagonal
            ija(k)=IA+1
          else
            sa(k)=e(i,j,2,1)  ! Store coupling on diagonal
            ija(k)=IA
          endif
!
        if(j.gt.2)then
          if(m.eq.1)then
            k=k+1
            sa(k)=d(i,j,1,1)  !  Store j-1 coefficient
            ija(k)=IA-2*JM2
            k=k+1
            sa(k)=d(i,j,1,2)  !  Store j-1 coefficient
            ija(k)=IA-2*JM2+1
          else
            k=k+1
            sa(k)=d(i,j,2,1)  !  Store j-1 coefficient
            ija(k)=IA-2*JM2
            k=k+1
            sa(k)=d(i,j,2,2)  !  Store j-1 coefficient
            ija(k)=IA-2*JM2+1
          endif
        endif
        if(i.gt.2)then
          if(m.eq.1)then
            k=k+1
            sa(k)=b(i,j,1,1)  !  Store i-1 coefficient
            ija(k)=IA-2
            k=k+1
            sa(k)=b(i,j,1,2)  !  Store i-1 coefficient
            ija(k)=IA-1
          else
            k=k+1
            sa(k)=b(i,j,2,1)  !  Store i-1 coefficient
            ija(k)=IA-2
            k=k+1
            sa(k)=b(i,j,2,2)  !  Store i-1 coefficient
            ija(k)=IA-1
          endif
        endif
        if(i.lt.JM1)then
          if(m.eq.1)then
            k=k+1
            sa(k)=a(i,j,1,1)  !  Store i+1 coefficient
            ija(k)=IA+2
            k=k+1
            sa(k)=a(i,j,1,2)  !  Store i+1 coefficient
            ija(k)=IA+3
          else
            k=k+1
            sa(k)=a(i,j,2,1)  !  Store i+1 coefficient
            ija(k)=IA+2
            k=k+1
            sa(k)=a(i,j,2,2)  !  Store i+1 coefficient
            ija(k)=IA+3
          endif
        endif
        if(j.lt.JM1)then
          if(m.eq.1)then
            k=k+1
            sa(k)=c(i,j,1,1)  !  Store j+1 coefficient
            ija(k)=IA+2*JM2
            k=k+1
            sa(k)=c(i,j,1,2)  !  Store j+1 coefficient
            ija(k)=IA+2*JM2+1
          else
            k=k+1
            sa(k)=c(i,j,2,1)  !  Store j+1 coefficient
            ija(k)=IA+2*JM2
            k=k+1
            sa(k)=c(i,j,2,2)  !  Store j+1 coefficient
            ija(k)=IA+2*JM2+1
          endif
        endif
          if(m.eq.1)then
            ija(IA+1)=k+1  !  As each row is completed, store index to next
          else
            ija(IA+2)=k+1  !  As each row is completed, store index to next
          endif
        enddo              !  End of loop on m
      enddo                !  End of loop on rows
!     Initialize for PBCG
      k =0
      do j=2,JM1
        do i=2,JM1
          do m=1,2
            k=k+1
            if(m.eq.1)then
              xA(k)=s(i,j)
            else
              xA(k)=w(i,j)
            endif
            bA(k)=f(i,j,m)
          enddo
        enddo
      enddo
      call linbcg(NP,bA,xA,ITOL,TOL,ITMAX,iter,err)
      write(*,'(/1x,a,e15.6)') 'Estimated error:',err
      write(*,'(/1x,a,i6)') 'Iterations needed:',iter
      call dsprsax(sa,ija,x,bcmp,NP)
C     this is a double precision version of sprsax

      j=2
      i=1
      do IA=1,NN-1,2
        i=i+1
        if(i.gt.JM1)then
          i=2
          j=j+1
        endif
        do m=1,2
          if(m.eq.1)then
            s(i,j)=xA(IA)
          else
            w(i,j)=xA(IA+1)
          endif
        enddo
      enddo
!       Store boundary values of vorticity into w(i,j) array
      do j=2,JM1
        w(1,j)=-2.d0*(s(2,j)-s(1,j))/del2-2.d0*vtbx(1,j)/del
     *         -alf*(vnbx(1,j+1)-vnbx(1,j-1))/(2.d0*del)
        w(JMAX,j)=-2.d0*(S(JM1,j)-s(JMAX,j))/del2+2.d0*vtbx(2,j)/del
     *            -alf*(vnbx(2,j+1)-vnbx(2,j-1))/(2.d0*del)
      enddo
      do i=2,JM1
        w(i,1)=-2.d0*alf2*(s(i,2)-s(i,1))/del2+2.d0*alf*vtby(i,1)/del
     *                   +(vnby(i+1,1)-vnby(i-1,1))/(2.d0*del)
        w(i,JMAX)=-2.d0*alf2*(s(i,JM1)-s(i,JMAX))/del2
     *           -2.d0*alf*vtby(i,2)/del
     *                   +(vnby(i+1,2)-vnby(i-1,2))/(2.d0*del)
         enddo
	w(1,1)=0.5d0*(w(1,2)+w(2,1))
	w(JMAX,1)=0.5d0*(w(JMAX,2)+w(JM1,1))
	w(JMAX,JMAX)=0.5d0*(w(JMAX,JM1)+w(JM1,JMAX))
	w(1,JMAX)=0.5d0*(w(1,JM1)+w(2,JMAX))
!      Calculate velocity
      do i=2,JM1
       do j=2,JM1
         u(i,j)=alf*(s(i,j+1)-s(i,j-1))/(2.d0*del)
         v(i,j)=-(s(i+1,j)-s(i-1,j))/(2.d0*del)
       enddo
      enddo
      do j=1,JMAX
        u(1,j)=vnbx(1,j)
        u(JMAX,j)=vnbx(2,j)
        v(1,j)=vtbx(1,j)
        v(JMAX,j)=vtbx(2,j)
      enddo
      do i=1,JMAX
        u(i,1)=vtby(i,1)
        u(i,JMAX)=vtby(i,2)
        v(i,1)=vnby(i,1)
        v(i,JMAX)=vnby(i,2)
      enddo
      do j=1,JMAX
        write(13,'(1x,100e12.3)')(u(i,j), i=1,JMAX)
        write(14,'(1x,100e12.3)')(v(i,j), i=1,JMAX)
        write(19,'(1x,100e12.3)')(s(i,j), i=1,JMAX)
        write(20,'(1x,100e12.3)')(w(i,j), i=1,JMAX)
      enddo

      do k=1,ija(NN+1)
        write(12,'(1x,I5,2x,e12.3)') ija(k),sa(k)
      enddo 
      END