double precision function multLat2d(diag,x1Diag,y1Diag,
     *     nX,nY, nMax, v,i,j)
      implicit none

      integer i,j, nX, nY, nMax
      integer iPlus1, iMinus1
      real*8 diag,x1Diag,y1Diag,v,sum

      dimension v(0:nMax+1, 0:nMax+1)
      

c *** Function multiplies banded penta diagonal matrix to formed by scalars 
c *** diag (on diagonal), 
c *** b (on subdiagonal) and c (on superdiag) and returns the 
c *** value of the i-th component of the product.      


      iPlus1 = mod(i, nX)+1
      iMinus1 = mod(nX+i-2, nX)+1

      sum = 0.d0

c     If statements taken out
c     for speed
c      if(j .GT. 0) then

      sum = sum + v(i,j-1)

c      endif
c      if(j .LT. ny+1) then

      sum = sum + v(i,j+1)

c      endif

      sum = y1Diag*sum + diag*v(i,j) 
     *     + x1Diag*(v(iMinus1,j)+v(iPlus1,j))

c      write (*,*) iMinus1, ",", i,
c     *     ",", iPlus1, "-->",
c     *     sum


c      if (dabs(sum) .GT. 1.d0) then
c         write (*,*) "Sum is ",sum
c         write (*,*) "(i,j) is (", i, ",",j,")"
c         write (*,*) "v Vals at mask are "
c         write (*,*) v(i,j-1), v(i,j+1), 
c     *        v(iMinus1, j), v(iPlus1, j),
c     *        v(i,j)
c         pause "Should be zero"
c      endif

c     Reminder to allocate arrays of proper size for boundary
c     conditions

      

      multLat2d = sum
      
      
      return
      end