c     Contains functions for computing 
c     linear and non-linear portions
c     of derivative
c
c     If you want to change the PDE, do
c     so here
c
c     Recall the following functison
c     from iniValFuncs.f
c
c     function getViscosity()
c     real*8
c
c     function getForcingX(x, y)
c     real*8
c
c     function getForcingY(x, y)
c     real*8
c
cccccccccccccccccccccccccccccccccccccccc
c     subroutine GetLinTerm
cccccccccccccccccccccccccccccccccccccccc
c     Gets the linear component
c     of the time derivative
c     of omega
c
c     Make absolutely sure
c     that boundaries are filled in
c     for psi and w before calling this!
c
c     Does not touch boundaries of derivative
c     Make sure to take this into account
c     
c     Note that some arguments are not touched
c     This is out of a preference to give
c     the function the exact same signature as
c     the non-linear function
      subroutine GetLinTerm(
     *     psi, w, x,y,
     *     nX, nY, nMax, 
     *     nMaxOut,
     *     termOut)

      implicit none
      integer nX, nY, nMax, nMaxOut
      integer i,j
      integer iPlus1, iMinus1

      real*8 psi, w, termOut
      real*8 x,y
      real*8 dPsiDx, dPsiDy
      real*8 dWdX, dWdY
      real*8 deltaX, deltaY
      real*8 diag,diagx,diagy
      real*8 multLat2d
      real*8 getViscosity

      dimension w(0:nMax+1,0:nMax+1)
      dimension psi(0:nMax+1,0:nMax+1)
      dimension termOut(0:nMaxOut+1, 0:nMaxOut+1)
      dimension x(0:nX+1), y(0:nY+1)

      deltaX = x(1)-x(0)
      deltaY = y(1)-y(0)

      diagx = 1.d0/(deltaX*deltaX)
      diagy = 1.d0/(deltaY*deltaY)

      diag = -2.d0*diagx-2.d0*diagy

      diagx = diagx*getViscosity()
      diagy = diagy*getViscosity()
      diag  = diag *getViscosity()

      do j=1,nY
         do i=1,nX
            termOut(i,j) = 
     *           multLat2d(diag,diagx,diagy,
     *           nX,nY, nMax, w,i,j)

         enddo
      enddo

      return
      end

cccccccccccccccccccccccccccccccccccccccc
c     subroutine GetNonLinTerm
cccccccccccccccccccccccccccccccccccccccc
c     Gets the non-linear component
c     of the time derivative
c     of omega
c
c     Make absolutely sure
c     that boundaries are filled in
c     for psi and w before calling this!
c
c     Does not touch boundaries of derivative
c     Make sure to take this into account
c     
      subroutine GetNonLinTerm(
     *     psi, w, x,y,
     *     nX, nY, nMax, 
     *     nMaxOut,
     *     termOut)

      implicit none
      integer nX, nY, nMax, nMaxOut
      integer i,j
      integer iPlus1, iMinus1

      real*8 psi, w, termOut
      real*8 x,y
      real*8 dPsiDx, dPsiDy
      real*8 dWdX, dWdY
      real*8 oneOver2dx
      real*8 oneOver2dy

      real*8 getForcingX,
     *     getForcingY

      dimension w(0:nMax+1,0:nMax+1)
      dimension psi(0:nMax+1,0:nMax+1)
      dimension termOut(0:nMaxOut+1, 0:nMaxOut+1)
      dimension x(0:nX+1), y(0:nY+1)

      oneOver2dx = x(1)-x(0)
      oneOver2dx = 5.d-1/oneOver2dx

      oneOver2dy = y(1)-y(0)
      oneOver2dy = 5.d-1/oneOver2dy

      do j=1,nY
         do i=1,nX

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



            dPsiDx = oneOver2dx*
     *           (
     *           psi(iPlus1,j)-psi(iMinus1,j)
     *           )
            dPsiDy = oneOver2dy*
     *           (
     *           psi(i,j+1)-psi(i,j-1)
     *           )
            dWdX = oneOver2dx*
     *           (
     *           w(iPlus1,j)-w(iMinus1,j)
     *           )
            dWdY = oneOver2dy*
     *           (
     *           w(i,j+1)-w(i,j-1)
     *           )

            termOut(i,j) = 
     *           dPsiDx*dWdY -
     *           dPsiDy*dWdX

ccccccccccccccccccccccccccccccccccccccccccccc
c     This is minus the first y derivative of
c     the x component of the forcing term
c     which is viscosity times the third
c     y derivative of the constant desired
c     uBar
            termOut(i,j) = termOut(i,j) +
     *           oneOver2dy*(
     *           getForcingX(x(i),y(j-1)) -
     *           getForcingX(x(i),y(j+1)) )

         enddo
      enddo
      

      return
      end