function bess2d(ix,iy,x,y,f,im,jm)
c
c=======================================================================
c                                                                    ===
c  This routine performs a 16-point cubic bessel interpolation at    ===
c  the grid locations X and Y from a regularly grided 2-D field F.   ===
c                                                                    ===
c     On Input:                                                      ===
c                                                                    ===
c     F        field to interpolate from (real array)                ===
c     IM, JM   first and second dimensions of F (integer)            ===
c     IX, IY   SW corner  grid position (grid units; integer)        ===
c     X, Y     position in the box (grid units; real)                ===
c     On Output:                                                     ===
c                                                                    ===
c     BESS2D    interpolated value (real)                            ===
c                                                                    ===
c     Calls:      BESS1D                                             ===
c                                                                    ===
c=======================================================================
c
c-----------------------------------------------------------------------
c  Define global data.
c-----------------------------------------------------------------------
c
#include <cdefs.h>
#include <pconst.h>
c
c-----------------------------------------------------------------------
c  Define local data.
c-----------------------------------------------------------------------
c
      integer i1,i2,i3,i4,im,ix,iy,j1,j2,j3,j4,jm
      FLOAT
     *      f1,f2,f3,f4,x,y
      FLOAT
     *      bess1d,bess2d
      FLOAT
     *      f(im,jm)
c
c-----------------------------------------------------------------------
c  Begin executable code.
c-----------------------------------------------------------------------
c
c-----------------------------------------------------------------------
c  Get the indices for the 4 x 4 raw data window.
c-----------------------------------------------------------------------
c
      i1=ix-1
      i2=ix
      i3=ix+1
      i4=ix+2
      j1=iy-1
      j2=iy
      j3=iy+1
      j4=iy+2
c
c-----------------------------------------------------------------------
c  Achieve the two-dimensional interpolation by a series of
c  one-dimensional Bessel interpolations.
c-----------------------------------------------------------------------
c
      f1=bess1d(y,f(i1,j1),f(i1,j2),f(i1,j3),f(i1,j4))
      f2=bess1d(y,f(i2,j1),f(i2,j2),f(i2,j3),f(i2,j4))
      f3=bess1d(y,f(i3,j1),f(i3,j2),f(i3,j3),f(i3,j4))
      f4=bess1d(y,f(i4,j1),f(i4,j2),f(i4,j3),f(i4,j4))
      bess2d=bess1d(x,f1,f2,f3,f4)
c
      return
      end