SUBROUTINE coef1
!     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
	INCLUDE 'include1.inc'
      do 12 i=2,JM1
        do 11 j=2,JM1
         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
!      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
11      continue
12    continue
!     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-vbx(1,j)/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+vbx(2,j)/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*vby(i,1)/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*vby(i,2)/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
      RETURN
	  END