function bes1d(x,f1,f2,f3,f4,ibry)
c
c=======================================================================
c                                                                    ===
c  This function performs a 1-d  cubic Bessel interpolation at       ===
c  the point X given the value of function at four consecutive       ===
c  points (F1,F2,F3,F4).   The value of the function at X lies       ===
c  between F2 and F3.                                                ===
c                                                                    ===
c     On Input:                                                      ===
c                                                                    ===
c     F1,F2,F3,F4  field to interpolate from (real)                  ===
c     X            position at which the field is interpolated (grid ===
c                  units; real)                                      ===
c     IBRY         pointer that indicates if either f2 or f3 is a    ===
c                  boundary point.   (integer)                       ===
c                    IBRY = 0  ->  neither f2 nor f3 is at a boundary===
c                    IBRY = 1  ->  f2 is at the boundary.            ===
c                    IBRY = 2  ->  f3 is at the boundary.            ===
c                                                                    ===
c     On Output:                                                     ===
c                                                                    ===
c     BES1D        interpolated value (real)                         ===
c                                                                    ===
c     Calls:          none                                           ===
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 ibry
      FLOAT bes1d
      FLOAT
     *      a0,a1,a2,a3,f1,f2,f3,f4,fp0,fp1,x
c
c-----------------------------------------------------------------------
c  Begin executable code.
c-----------------------------------------------------------------------
c
      if((ibry.lt.1).or.(ibry.gt.2)) then
        fp0=p5*(f3-f1)
        fp1=p5*(f4-f2)
      elseif(ibry.eq.1) then
        fp0=p5*(-c3*f2+c4*f3-f4)
        fp1=p5*(f4-f2)
      else
        fp0=p5*(f3-f1)
        fp1=p5*(c3*f3-c4*f2+f1)
      endif
      a0=f2
      a1=fp0
      a2=c3*f3-fp1-c3*a0-c2*a1
      a3=f3-a0-a1-a2
c
c-----------------------------------------------------------------------
c  Interpolate F(X) with a cubic polynomial.
c-----------------------------------------------------------------------
c
      bes1d=a0+x*(a1+x*(a2+a3*x))
      return
      end