SUBROUTINE BC3
	INCLUDE 'include3.inc'
!      Aspect ratio and boundary conditions, velocities either 0 or 1
!      icase=1: Cartesian coordinates
!       ibc=0    Zero normal components of velocity
!       ibc=1    Couette flow
!       ibc=2    Plane Poiseuille
!       ibc=3    Boundary layer flow along y=0 surface
!      icase=2: Cylindrical polar coordinates
!       ibc=1    Radial flow
!       ibc=2    Couette flow
!       ibc=3    Flow around a cylinder
!      nframes = 0 or up to 20 frames when streamlines are captured
      double precision sum
      read(8,*) icase,ibc,alf,beta,Re,tmax,dtmax,irestart
      read(8,*) vtbx1,vtbx2,vtby1,vtby2,vnbx1,vnbx2,vnby1,vnby2
	read(8,*) nframes
      if(nframes.gt.0)then
	  nframe=1
        read(8,*) (tim(n),n=1,nframes)
	  write(9,'(1x,50e12.4)') Re,(tim(n),n=1,nframes)
	endif
!      Parameter specification
	dwmax=5.d-1
	dtinit=1.d-4
	t=0.d0
	eps=1.d-3
	gamma=1.d0/dlog(beta)
	gam2=gamma**2
      alf2=alf**2
      del=1.d0/(JMAX-1)
      del2=del**2
	if(icase.eq.2)then
	  do i=1,JMAX
	    r(i)=dexp(del*(i-1)/gamma)
	  enddo
	endif
!     Initialize solution
      do i=1,JMAX
        do j=1,JMAX
          s(i,j)=0.d0
          w(i,j)=0.d0
          dw(i,j)=0.d0
	    p(i,j)=0.d0
        enddo
      enddo
	if(irestart.eq.1)then		! Restart from previous run
	do j=1,JMAX
        read(19,'(1x,100e12.3)')(s(i,j), i=1,JMAX)
        read(20,'(1x,100e12.3)')(w(i,j), i=1,JMAX)
	enddo
	REWIND 19
	REWIND 20
	endif
!     Boundary conditions
!     sbx & sby are boundary values of stream function
!     vtbx & vtby are boundary conditions on tangental components of VELOCITY
!     vnbx & vnby are boundary conditions on normal components of VELOCITY
      if(icase.eq.1)then    ! Cartesian coordinates
!     x boundaries, i.e., left and right, x=0,1
      do j=1,JMAX
	  if(ibc.eq.0)then
	    vnbx(1,j)=vnbx1                                   ! Constant
	    vnbx(2,j)=vnbx2
	  elseif(ibc.eq.1)then
          vnbx(1,j)=vnbx1*del*(j-1)	                      ! Linear profile
          vnbx(2,j)=vnbx2*del*(j-1)	                      ! Linear profile
	  elseif(ibc.eq.2)then
          vnbx(1,j)=4.d0*vnbx1*(del*(j-1)-(del*(j-1))**2) ! Parabolic profile
          vnbx(2,j)=4.d0*vnbx2*(del*(j-1)-(del*(j-1))**2) ! Parabolic profile
	  elseif(ibc.eq.3)then								! Boundary layer
!	    if(j.eq.1)then
!	       vnbx(1,j)=0.5d0
!	    else
 		   vnbx(1,j)=1.d0
!	    endif                                  
 	    vnbx(2,j)=erf(((j-1)*del/alf)/sqrt(8.d0/Re))    ! Boundary layer
	  endif
       vtbx(1,j)=vtbx1
       vtbx(2,j)=vtbx2
      enddo
	if(ibc.eq.3)then                                    ! Boundary layer
	sum=0.d0
	do j=2,JMAX
	  sum=sum+0.5d0*del*(vnbx(1,j-1)+vnbx(1,j))/alf		! Change in flux on x boundaries
     *         -0.5d0*del*(vnbx(2,j-1)+vnbx(2,j))/alf
	enddo
	do j=2,JMAX											! Transverse velocity
	  vtbx(1,j)=0.d0
	  vtbx(2,j)=2.d0*sum*(j-1)/JM1
	enddo
	endif
!     y boundaries,i.e., y=0,1
      do i=1,JMAX
       vtby(i,1)=vtby1
       vtby(i,2)=vtby2
	  if(ibc.eq.0)then
	    vnby(i,1)=vnby1                                   !  Constant
	    vnby(i,2)=vnby2
	  elseif(ibc.eq.1)then
          vnby(i,1)=vnby1*del*(i-1)	                      ! Linear profile
          vnby(i,2)=vnby2*del*(i-1)	                      ! Linear profile
	  elseif(ibc.eq.2)then
          vnby(i,1)=4.d0*vnby1*(del*(i-1)-(del*(i-1))**2)	  !  Parabolic profile
          vnby(i,2)=4.d0*vnby2*(del*(i-1)-(del*(i-1))**2)	  !  Parabolic profile
	  elseif(ibc.eq.3)then								  ! Boundary layer
	    if(i.eq.1)then
	      vtby(i,1)=0.0d0
	  else
	      vtby(i,1)=0.d0
	  endif
	    vnby(i,1)=0.d0                                    
          vnby(i,2)=2.d0*(i-1)*sum/JM1
	  endif
      enddo
!     Integrate clockwise along boundaries to calculate B.C. for stream function
      sbx(1,1)=0.d0
	do j=2,JMAX
	  sbx(1,j)=sbx(1,j-1)+0.5d0*del*(vnbx(1,j-1)+vnbx(1,j))/alf
	enddo
	sby(1,2)=sbx(1,JMAX)
	do i=2,JMAX
	  sby(i,2)=sby(i-1,2)-0.5d0*del*(vnby(i-1,2)+vnby(i,2))
	enddo
	sbx(2,JMAX)=sby(JMAX,2)
	do j=JM1,1,-1
	  sbx(2,j)=sbx(2,j+1)-0.5d0*del*(vnbx(2,j+1)+vnbx(2,j))/alf
	enddo
	sby(JMAX,1)=sbx(2,1)
	do i=JM1,1,-1
	  sby(i,1)=sby(i+1,1)+0.5d0*del*(vnby(i+1,1)+vnby(i,1))
	enddo
	else            ! Cylindrical polar coordinates
!     r boundaries, i.e.,  r=1,beta; z=0,1
      do j=1,JMAX
	  if(ibc.eq.1)then
          vnbx(1,j)=vnbx1/r(1)    	                      ! Radial flow
          vnbx(2,j)=vnbx2/r(JMAX)	                          
          vtbx(1,j)=vtbx1
          vtbx(2,j)=vtbx2
	  elseif(ibc.eq.2)then
          vnbx(1,j)=0.d0                                    ! Couette flow
          vnbx(2,j)=0.d0									  
          vtbx(1,j)=vtbx1
          vtbx(2,j)=vtbx2
	  else										! potential flow past cylinder
	    vnbx(1,j)=0.d0                                    	! no flux
		vnbx(2,j)=+dcos(del*(j-1)/alf)*(1.d0-r(JMAX)**(-2))	! Potential flow
		vtbx(1,j)=0.d0				                        ! no slip	!!!!
		vtbx(2,j)=-dsin(del*(j-1)/alf)*(1.d0+r(JMAX)**(-2))	! Potential flow                                    
	  endif
      enddo
!     theta boundaries,i.e., theta=0,1
      do i=1,JMAX
	  if(ibc.eq.1)then
          vnby(i,1)=0.d0        	                          ! Radial flow
          vnby(i,2)=0.d0      	                          
          vtby(i,1)=vtby1/r(i)
          vtby(i,2)=vtby2/r(i)
	  elseif(ibc.eq.2)then
          vnby(i,1)=(r(i)-1.d0/r(i))/(beta-1.d0/beta)       ! Couette flow
          vnby(i,2)=(r(i)-1.d0/r(i))/(beta-1.d0/beta)       
	    vtby(i,1)=0.d0
	    vtby(i,2)=0.d0
	  else												  ! flow past cylinder
	    vnby(i,1)=0.d0                 ! Symmetry                   
		vnby(i,2)=0.d0
		vtby(i,1)=+(1.d0-r(i)**(-2))	   ! BC based on potential flow
		vtby(i,2)=+dcos(1.d0/alf)*(1.d0-r(i)**(-2)) 
	  endif
      enddo
!     Integrate clockwise along boundaries to calculate B.C. for stream function
      sbx(1,1)=0.d0		   ! r = 1; theta = 0
	do j=2,JMAX			   ! r = 1 boundary
	  sbx(1,j)=sbx(1,j-1)+0.5d0*r(1)*del*(vnbx(1,j-1)+vnbx(1,j))/alf
	enddo
	sby(1,2)=sbx(1,JMAX)
	do i=2,JMAX			   ! theta = pi  boundary
	  sby(i,2)=sby(i-1,2)
     *  -0.5d0*(vnby(i-1,2)+vnby(i,2))*(r(i)-r(i-1))
	enddo
	sbx(2,JMAX)=sby(JMAX,2)
	do j=JM1,1,-1		   ! r = beta boundary
	 sbx(2,j)=sbx(2,j+1)-0.5d0*r(JMAX)*del*(vnbx(2,j+1)+vnbx(2,j))/alf
	enddo
	sby(JMAX,1)=sbx(2,1)
	do i=JM1,1,-1		   ! theta = 0 boundary
	  sby(i,1)=sby(i+1,1)
     *  -0.5d0*(vnby(i+1,1)+vnby(i,1))*(r(i)-r(i+1))
	enddo
	endif
!	Following for either coordinates
!     Load boundary values of stream function into s(i,j) array
      do j=1,JMAX
        s(1,j)=sbx(1,j)
        s(JMAX,j)=sbx(2,j)
      enddo
      do i=1,JMAX
        s(i,1)=sby(i,1)
        s(i,JMAX)=sby(i,2)
      enddo
      RETURN
	END