SUBROUTINE coef3
	INCLUDE 'include3.inc'
!
!     Calculate coefficients, e.g.,
!     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
!     Index 1 for stream function; index 2 for vorticity, diagonal terms
      do i=2,JM1
        do j=2,JM1
	   if(icase.eq.1)then  ! Cartesian coordinates
           a(i,j,1,1)=1.d0
           b(i,j,1,1)=1.d0
           c(i,j,1,1)=alf2
           d(i,j,1,1)=alf2
           e(i,j,1,1)=-2.d0*(1.d0+alf2)
           f(i,j,1)=0.d0
           a(i,j,2,2)=1.d0
           b(i,j,2,2)=1.d0
           c(i,j,2,2)=alf2
           d(i,j,2,2)=alf2
           e(i,j,2,2)=-2.d0*(1.d0+alf2)
           f(i,j,2)=0.d0
	   else				   ! Cylindrical polar coordinates
           a(i,j,1,1)=gam2/(r(i)**2)
           b(i,j,1,1)=gam2/(r(i)**2)
           c(i,j,1,1)=alf2/(r(i)**2)
           d(i,j,1,1)=alf2/(r(i)**2)
           e(i,j,1,1)=-(a(i,j,1,1)+b(i,j,1,1)+c(i,j,1,1)+d(i,j,1,1))
           f(i,j,1)=0.d0
           a(i,j,2,2)=gam2/(r(i)**2)
           b(i,j,2,2)=gam2/(r(i)**2)
           c(i,j,2,2)=alf2/(r(i)**2)
           d(i,j,2,2)=alf2/(r(i)**2)
           e(i,j,2,2)=-(a(i,j,1,1)+b(i,j,1,1)+c(i,j,1,1)+d(i,j,1,1))
           f(i,j,2)=0.d0
	   endif
!      Coupling of dependent variables between equations
         a(i,j,1,2)=0.d0
         b(i,j,1,2)=0.d0
         c(i,j,1,2)=0.d0
         d(i,j,1,2)=0.d0
         e(i,j,1,2)=del2    !  Coupling of stream function with vorticity
         a(i,j,2,1)=0.d0
         b(i,j,2,1)=0.d0
         c(i,j,2,1)=0.d0
         d(i,j,2,1)=0.d0
         e(i,j,2,1)=0.d0
        enddo
	enddo
!     Calculate CONVECTIVE TERMS, look for upstream direction
      do j=2,JM1
        do i=2,JM1
	  if(icase.eq.1)then			   ! Cartesian coordinates
          uim1=0.5d0*(u(i,j)+u(i-1,j))   ! u at i-1/2
	    if(uim1.gt.0.)then
	      b(i,j,2,2)=b(i,j,2,2)+del*Re*uim1
	    else
	      e(i,j,2,2)=e(i,j,2,2)+del*Re*uim1
	    endif
          uip1=0.5d0*(u(i,j)+u(i+1,j))   ! u at i+1/2
          if(uip1.gt.0.)then
	      e(i,j,2,2)=e(i,j,2,2)-del*Re*uip1
	    else
	      a(i,j,2,2)=a(i,j,2,2)-del*Re*uip1
	    endif
	    vjm1=0.5d0*(v(i,j)+v(i,j-1))    ! v at j-1/2
	    if(vjm1.gt.0.)then
	      d(i,j,2,2)=d(i,j,2,2)+alf*del*Re*vjm1
	    else
	      e(i,j,2,2)=e(i,j,2,2)+alf*del*Re*vjm1
	    endif
	    vjp1=0.5d0*(v(i,j)+v(i,j+1))   ! v at j+1/2
	    if(vjp1.gt.0.)then
	      e(i,j,2,2)=e(i,j,2,2)-alf*del*Re*vjp1
	    else
	      c(i,j,2,2)=c(i,j,2,2)-alf*del*Re*vjp1
	    endif
	  else							  ! cylindrical polar
          uim1=0.5d0*(u(i,j)+u(i-1,j))   ! u at i-1/2
	    if(uim1.gt.0.)then
	      b(i,j,2,2)=b(i,j,2,2)
     *	          +gamma*del*Re*uim1*(r(i)+r(i-1))/(2.d0*r(i)**2)
	    else
	      e(i,j,2,2)=e(i,j,2,2)
	*              +gamma*del*Re*uim1*(r(i)+r(i-1))/(2.d0*r(i)**2)
	    endif
          uip1=0.5d0*(u(i,j)+u(i+1,j))   ! u at i+1/2
          if(uip1.gt.0.)then
	      e(i,j,2,2)=e(i,j,2,2)
     *              -gamma*del*Re*uip1*(r(i)+r(i+1))/(2.d0*r(i)**2)
	    else
	      a(i,j,2,2)=a(i,j,2,2)
	*              -gamma*del*Re*uip1*(r(i)+r(i+1))/(2.d0*r(i)**2)
	    endif
	    vjm1=0.5d0*(v(i,j)+v(i,j-1))    ! v at j-1/2
	    if(vjm1.gt.0.)then
	      d(i,j,2,2)=d(i,j,2,2)+alf*del*Re*vjm1/r(i)
	    else
	      e(i,j,2,2)=e(i,j,2,2)+alf*del*Re*vjm1/r(i)
	    endif
	    vjp1=0.5d0*(v(i,j)+v(i,j+1))   ! v at j+1/2
	    if(vjp1.gt.0.)then
	      e(i,j,2,2)=e(i,j,2,2)-alf*del*Re*vjp1/r(i)
	    else
	      c(i,j,2,2)=c(i,j,2,2)-alf*del*Re*vjp1/r(i)
	    endif
	  endif
	  enddo
	enddo
!	BOUNDARY CONDITIONS
      if(icase.eq.1)then         ! Cartesian coordinates 
!     x boundaries,i.e. left and right

      do i=2,JM1,JM3
        do j=2,JM1
          if(i.eq.2)then
            f(i,j,1)=f(i,j,1)-b(i,j,1,1)*sbx(1,j)
            f(i,j,2)=f(i,j,2)
     *                 -b(i,j,2,2)*(2.d0*(sbx(1,j)/del2-vtbx(1,j)/del)
     *                 -alf*(vnbx(1,j+1)-vnbx(1,j-1))/(2.d0*del))
            e(i,j,2,1)=e(i,j,2,1)-b(i,j,2,2)*2.d0/del2
            b(i,j,1,1)=0.d0
            b(i,j,2,2)=0.d0
          endif
          if(i.eq.JM1)then
            f(i,j,1)=f(i,j,1)-a(i,j,1,1)*sbx(2,j)
            f(i,j,2)=f(i,j,2)
     *                 -a(i,j,2,2)*(2.d0*(sbx(2,j)/del2+vtbx(2,j)/del)
     *                 -alf*(vnbx(2,j+1)-vnbx(2,j-1))/(2.d0*del))
            e(i,j,2,1)=e(i,j,2,1)-a(i,j,2,2)*2.d0/del2
            a(i,j,1,1)=0.d0
            a(i,j,2,2)=0.d0
          endif
        enddo
      enddo

!     y boundaries, i.e., top and bottom

      do j=2,JM1,JM3
        do i=2,JM1
          if(j.eq.2)then
            f(i,j,1)=f(i,j,1)-d(i,j,1,1)*sby(i,1)
            f(i,j,2)=f(i,j,2)
     *          -d(i,j,2,2)*(2.d0*(alf2*sby(i,1)/del2+alf*vtby(i,1)/del)
     *                   +(vnby(i+1,1)-vnby(i-1,1))/(2.d0*del))  
            e(i,j,2,1)=e(i,j,2,1)-d(i,j,2,2)*2.d0*alf2/del2
            d(i,j,1,1)=0.d0
            d(i,j,2,2)=0.d0
         endif
         if(j.eq.JM1)then
           f(i,j,1)=f(i,j,1)-c(i,j,1,1)*sby(i,2)
           f(i,j,2)=f(i,j,2)
     *          -c(i,j,2,2)*(2.d0*(alf2*sby(i,2)/del2-alf*vtby(i,2)/del)
     *                   +(vnby(i+1,2)-vnby(i-1,2))/(2.d0*del))  
           e(i,j,2,1)=e(i,j,2,1)-c(i,j,2,2)*2.d0*alf2/del2
           c(i,j,1,1)=0.d0
           c(i,j,2,2)=0.d0
         endif
        enddo
      enddo
	else               !   Cylindrical polar coordinates
!     r boundaries,i.e. inner and outer radaii

      do i=2,JM1,JM3
        do j=2,JM1
          if(i.eq.2)then
            f(i,j,1)=f(i,j,1)-b(i,j,1,1)*sbx(1,j)
            f(i,j,2)=f(i,j,2)
     *        -b(i,j,2,2)*(2.d0*(gam2*sbx(1,j)/(del2*r(1)**2)				  	  
     *                          -gamma*vtbx(1,j)/(r(1)*del))
     *                 -alf*(vnbx(1,j+1)-vnbx(1,j-1))/(2.d0*r(1)*del))
            e(i,j,2,1)=e(i,j,2,1)-b(i,j,2,2)*2.d0*gam2/(del2*r(1)**2)			  
            b(i,j,1,1)=0.d0
            b(i,j,2,2)=0.d0
          endif
          if(i.eq.JM1)then
             f(i,j,1)=f(i,j,1)-a(i,j,1,1)*sbx(2,j)
             a(i,j,1,1)=0.d0
             f(i,j,2)=f(i,j,2)
     *        -a(i,j,2,2)*(2.d0*(gam2*sbx(2,j)/(del2*r(JMAX)**2)				  
     *                          +gamma*vtbx(2,j)/(r(JMAX)*del))
     *                -alf*(vnbx(2,j+1)-vnbx(2,j-1))/(2.d0*r(JMAX)*del))
            e(i,j,2,1)=e(i,j,2,1)-a(i,j,2,2)*2.d0*gam2/(del2*r(JMAX)**2)		  
            a(i,j,2,2)=0.d0
          endif
        enddo
      enddo

!     theta boundaries, i.e., theta = 0, pi
	if(ibc.lt.3)then     ! radial or Couette flow
      do j=2,JM1,JM3
        do i=2,JM1
          if(j.eq.2)then
            f(i,j,1)=f(i,j,1)-d(i,j,1,1)*sby(i,1)
            f(i,j,2)=f(i,j,2)
     *          -d(i,j,2,2)*(2.d0*(alf2*sby(i,1)/(del2*r(i)**2)	   
     *                            +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))  
            e(i,j,2,1)=e(i,j,2,1)-d(i,j,2,2)*2.d0*alf2/(del2*r(i)**2)
            d(i,j,1,1)=0.d0
            d(i,j,2,2)=0.d0
         endif
         if(j.eq.JM1)then
           f(i,j,1)=f(i,j,1)-c(i,j,1,1)*sby(i,2)
           f(i,j,2)=f(i,j,2)
     *          -c(i,j,2,2)*(2.d0*(alf2*sby(i,2)/(del2*r(i)**2)	   
     *                            -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))  
           e(i,j,2,1)=e(i,j,2,1)-c(i,j,2,2)*2.d0*alf2/(del2*r(i)**2)
           c(i,j,1,1)=0.d0
           c(i,j,2,2)=0.d0
         endif
        enddo
	enddo
	else              ! flow past cylinder, zero vorticity at theta = 0, pi
      do j=2,JM1,JM3
        do i=2,JM1
          if(j.eq.2)then
            f(i,j,1)=f(i,j,1)-d(i,j,1,1)*sby(i,1)
            d(i,j,1,1)=0.d0
            d(i,j,2,2)=0.d0
         endif
         if(j.eq.JM1)then
           f(i,j,1)=f(i,j,1)-c(i,j,1,1)*sby(i,2)
           c(i,j,1,1)=0.d0
           c(i,j,2,2)=0.d0
         endif
        enddo
      enddo
	endif  ! end of ibc
	endif  ! end of icase
!
!     Update coefficients for transient solution
      do i=2,JM1
	  do j=2,JM1
          f(i,j,1)=f(i,j,1)-e(i,j,1,2)*w(i,j)
	    f(i,j,2)=f(i,j,2)-a(i,j,2,2)*w(i+1,j)-b(i,j,2,2)*w(i-1,j)
     *     -c(i,j,2,2)*w(i,j+1)-d(i,j,2,2)*w(i,j-1)-e(i,j,2,2)*w(i,j)
          e(i,j,2,2)=e(i,j,2,2)-del2*Re/dt
        enddo
      enddo
      RETURN
	END