PROGRAM STREAMF
! Stream function and vorticity solution of Navier-Stokes equation in 2-D
! Transient and finite Reynolds number 	
! s(i,j) is the z component of the vector potential or the stream function
! w(i,j) is the z component of the vorticity
! dw(i,j) is the change in w over the time step
! u(i,j) and v(i,j) are x and y components of velocity, respectively
! (x,y)=(xi,yj) or r and theta components in cylindrical polar coordinates
! 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
! p is pressure and tau is shear stress
! icase: 1 = Cartesian coordinates (x,y); 
!        2 = cylindrical polar [z=ln(r)/ln(r2/r1),theta]
! ibc: For Cartesian coordinates 1 linear(Couette); 2=parabolic(plane-Poiseuille); 3 constant
! ibc: For cylindrical polar,1=radial flow;2=Couette flow, 3=flow past cylinder
! beta = r2/r1, gamma=1/ln(beta)
! alf = 1/pi or 1/reference angle in cylindrical coordinates; 0<theta<1
! r is dimensionless radius, 1<r<beta
! z=gamma*ln(r); 0<z<1
!
! 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
! eps is epsilon for pressure equation
! PBCG from Numerical Recipes is used to solve linear system of equations

	INCLUDE 'include3.inc'
	OPEN( 8, FILE='input.dat')
	OPEN( 9, FILE='pram.dat')
      OPEN(12, FILE='matrix.dat')
      OPEN(13, FILE='u.dat')
      OPEN(14, FILE='v.dat')
      OPEN(15, FILE='p.dat')
	OPEN(16, FILE='f.dat')
	OPEN(17, FILE='tau.dat')
      OPEN(19, FILE='s.dat')
      OPEN(20, FILE='w.dat')
	OPEN(21, FILE='s1.dat')
	OPEN(22, FILE='s2.dat')
	OPEN(23, FILE='s3.dat')
	OPEN(24, FILE='s4.dat')
	OPEN(25, FILE='s5.dat')
	OPEN(26, FILE='s6.dat')
	OPEN(27, FILE='s7.dat')
	OPEN(28, FILE='s8.dat')
	OPEN(29, FILE='s9.dat')
      if(NMAX.lt.10*JMAX2)then
        write(*,'(/1x,a)') 'NMAX is too small. Redimension'
        stop
      endif
      CALL bc3
!
!     Start of time loop
      dt=dtinit
	do 100 while (t.lt.tmax)
	t=t+dt
	  do 100 it=1,2
	CALL coef3
	CALL store3
!     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)=dw(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
	dwm=0.d0
      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
            dw(i,j)=xA(IA+1)
	      if(it.eq.2)then          ! Only update after corrector step
              w(i,j)=w(i,j)+dw(i,j)  ! Update vorticity
              if(abs(dw(i,j)).gt.dwm) dwm=abs(dw(i,j))  ! Maximum dw
            endif
          endif
        enddo
      enddo
!      Calculate velocity
      do i=2,JM1
       do j=2,JM1
	   if(icase.eq.1)then	   ! Cartesian coordinates
           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)
	   else                    ! Cylindrical polar coordinates
	     u(i,j)=alf*(s(i,j+1)-s(i,j-1))/(2.d0*r(i)*del)
	     v(i,j)=-gamma*(s(i+1,j)-s(i-1,j))/(2.d0*r(i)*del)
	   endif
       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

!       Time step control
      if(it.eq.2)then
	 if(nframes.eq.0)then
	  if(t.gt.0.99*tmax.and.t.le.tmax) then
	    t=1.00001d0*tmax
	  else
	    dt=dt*dwmax/(dwm+1.d-6)
	    if(dt.gt.dtmax) dt=dtmax   ! Limit time step size to see oscillations
	    if(t+dt.gt.tmax) dt=0.99999999*(tmax-t)
	  endif
	 else        !  Capture frames for movie
	  if(t.gt.0.99*tim(nframe).and.t.lt.tim(nframe)) then
	    t=1.000001*tim(nframe)
	    nfile=20+nframe
	    do j=1,JMAX
	    write(nfile,'(1x,100e12.3)')(s(i,j), i=1,JMAX)
	    enddo
	    nframe=nframe+1
	    if(nframe.gt.nframes) nframes=0
	    dt=dtold
	    if(nframes.ne.0.AND.t+dt.gt.tim(nframe))
     *		 dt=0.9999999*(tim(nframe)-t)
	  else
	    dt=dt*dwmax/dwm
	    if(dt.gt.dtmax) dt=dtmax   ! Limit time step size to see oscillations
	    dtold=dt
	    if(t+dt.gt.tim(nframe)) dt=0.9999999*(tim(nframe)-t)
	  endif
	 endif
	  write(*,'(/1x,a,3e10.2)') 'Time, dt, dwm = ', t, dt, dwm
	endif
100   enddo	! End of time loop
!       Store boundary values of vorticity into w(i,j) array
      if(icase.eq.1)then             ! Cartesian coordinates
      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
	else                          ! Cylindrical polar coordinates
      do j=2,JM1					  ! Radial boundaries
        w(1,j)=-2.d0*gam2*(s(2,j)-s(1,j))/(del2*r(1)**2)
     *         -2.d0*gamma*vtbx(1,j)/(del*r(1))
     *         -alf*(vnbx(1,j+1)-vnbx(1,j-1))/(2.d0*del*r(1))
        w(JMAX,j)=-2.d0*gam2*(s(JM1,j)-s(JMAX,j))/(del2*r(JMAX)**2)
     *            +2.d0*gamma*vtbx(2,j)/(del*r(JMAX))
     *            -alf*(vnbx(2,j+1)-vnbx(2,j-1))/(2.d0*r(JMAX)*del)
      enddo
      do i=2,JM1					  ! Theta boundaries
	    if(ibc.lt.3)then
        w(i,1)=-2.d0*alf2*(s(i,2)-s(i,1))/(del2*r(i)**2)
     *	  +2.d0*alf*vtby(i,1)/(r(i)*del)
     * +gamma*(r(i+1)*vnby(i+1,1)-r(i-1)*vnby(i-1,1))/(2.d0*del*r(i)**2)  
        w(i,JMAX)=-2.d0*alf2*(s(i,JM1)-s(i,JMAX))/(del2*r(i)**2)
     *      -2.d0*alf*vtby(i,2)/(r(i)*del)
     * +gamma*(r(i+1)*vnby(i+1,2)-r(i-1)*vnby(i-1,2))/(2.d0*del*r(i)**2)
	    else  !Flow past cylinder for theta = 0&pi has w=0 because of symmetry
		w(i,1)=0.d0
		w(i,JMAX)=0.d0
		endif  
      enddo
	endif
	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))
      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
!     Calculate pressure distribution
	CALL press3
      do j=1,JMAX
        write(15,'(1x,100e12.4)')(p(i,j), i=1,JMAX)
        write(16,'(1x,100e12.4)')(f(i,j,1), i=1,JMAX)
        write(17,'(1x,100e12.4)')(tau(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