double precision function approxit(f,mmax,nmax, m, n, a, b, c, d, x ,y)

      implicit none

      integer m,n,mmax,nmax
      double precision x,y,a,b,c,d,f(mmax+2,nmax+2)

      double precision u,v
      integer i,j

      i = int((m+1)*(x-a)/(b-a)+1)
      j = int((n+1)*(y-c)/(d-c)+1)

      if(i .lt. 1) then
c         write(*,*)'Out of range, x < a.'
         i = 1
      else if(i .gt. m+1) then
c         write(*,*)'Out of range, x > b.'
         i = m+1
      endif      
      if(j .lt. 1) then
c         write(*,*)'Out of range, y < c.'
         j = 1
      else if(j .gt. n+1) then
c         write(*,*)'Out of range, y > d.'
         j = n+1
      endif
      u = (m+1)*(x-a)/(b-a) - i +1
      v = (n+1)*(y-c)/(d-c) - j +1
      approxit = (f(i,j)*u+f(i+1,j)*(1-u))*v+(f(i,j+1)*u+f(i+1,j+1)*(1-u))*(1-v)
      
      return
      end