SUBROUTINE bc2
	INCLUDE 'include2.inc'
!       ibc=1 for Coulette flow; ibc=2 for Hele-Shaw
	read(8,*) ibc,alf,vtbx1,vtbx2,vtby1,vtby2,vnbx1,vnbx2,vnby1,vnby2
      alf2=alf**2
      del=1.d0/(JMAX-1)
      del2=del**2
!     Initialize solution
      do i=1,JMAX
        do j=1,JMAX
          s(i,j)=0.d0
          w(i,j)=0.d0
        enddo
      enddo
!     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
!     x boundaries, i.e., left and right, x=0,1
      do j=1,JMAX
	if(ibc.eq.1)then
       vnbx(1,j)=vnbx1*del*(j-1)	                      ! Linear profile
       vnbx(2,j)=vnbx2*del*(j-1)	                      ! Linear profile
      else
        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
	endif
       vtbx(1,j)=vtbx1
       vtbx(2,j)=vtbx2
      enddo
!     y boundaries,i.e., y=0,1
      do i=1,JMAX
	if(ibc.eq.1)then
       vnby(i,1)=vnby1*del*(i-1)	                          ! Linear profile
       vnby(i,2)=vnby2*del*(i-1)	                          ! Linear profile
	else
       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
	endif
       vtby(i,1)=vtby1
       vtby(i,2)=vtby2
      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.5*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.5*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.5*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.5*del*(vnby(i+1,1)+vnby(i,1))
	enddo
!     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