subroutine bilinear(x,y,x1,x2,y1,y2,f1,f2,f3,f4,fun)

      real x,y,f1,f2,f3,f4,fun
      real x1,x2,y1,y2,t,u

* *** ********** Bilinear interpolation of 2D grid ******************
* *** ***************************************************************
*     fun = f(x,y)
*     (x,y) are within two-dimensional grid :

*           (4) f(J,K+1)      (3) f(J+1,K+1)    y2
*                       (x,y)
*           (1) f(J,K)        (2) f(J+1,K)      y1
*
*                 x1               x2
* *** ***************************************************************

      t = (x-x1)/(x2-x1)
      u = (y-y1)/(y2-y1)

      fun = (1.-t)*(1.-u)*f1  + t*(1.-u)*f2 + t*u*f3 + (1.-t)*u*f4

      end

* *******************************************************************
* ***

      subroutine linear(x,x1,x2,f1,f2,f)
      real x,x1,x2,f1,f2,f

      t= (x-x1)/(x2-x1)

      f = f1 + t*(f2-f1)

      end

* *** **************************************************************
* ***