function bess2d_msk (ix,iy,x,y,f,mask,im,jm,act_val,flag_val)
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  or approximations thereof for a masked field.                     ===
c                                                                    ===
c  ------                                                            ===
c  Input:                                                            ===
c  ------                                                            ===
c                                                                    ===
c     ACT_VAL      active value of mask.        (integer)            ===
c     F            field to interpolate from    (real array)         ===
c     IM, JM       1st & 2nd dimensions of F    (integer)            ===
c     MASK         mask array.                  (integer array)      ===
c     IX, IY       SW corner  grid position     (integer; grid units)===
c     X, Y         position in the box          (real; grid units)   ===
c                                                                    ===
c  -------                                                           ===
c  Output:                                                           ===
c  -------                                                           ===
c                                                                    ===
c     BESS2D_MSK   interpolated value           (real)               ===
c                                                                    ===
c  Calls:  BESS2D,  EXTRAP2                                          ===
c                                                                    ===
c=======================================================================
c
c-----------------------------------------------------------------------
c  Define global data.
c-----------------------------------------------------------------------
c
#include <cdefs.h>
#include <pconst.h>
#include <iounits.h>
c
c-----------------------------------------------------------------------
c  Define local data.
c-----------------------------------------------------------------------
c
      integer act_val,i,im,ix,ixwk,iy,iywk,j,jm,k
      integer mask(im,jm),ms(4,4)
      logical cell_ok
      FLOAT
     *      flag_val,s,t,x,y
      FLOAT
     *      fs(4,4)
      FLOAT
     *      bess2d,bess2d_msk
      FLOAT
     *      f(im,jm)
      integer i0(3,8),j0(3,8),i1(3,4),j1(3,4)
c
      save i0,j0,i1,j1
c
      data i0/
     *     1,2,3,
     *     1,2,3,
     *     2,2,2,
     *     3,3,3,
     *     4,3,2,
     *     4,3,2,
     *     3,3,3,
     *     2,2,2/
      data j0/
     *     3,3,3,
     *     2,2,2,
     *     1,2,3,
     *     1,2,3,
     *     2,2,2,
     *     3,3,3,
     *     4,3,2,
     *     4,3,2/
      data i1/
     *     1,1,2,
     *     1,1,2,
     *     4,3,4,
     *     4,3,4/
      data j1/
     *     4,3,4,
     *     1,2,1,
     *     1,1,2,
     *     4,4,3/
c
c=======================================================================
c  Begin executable code.
c=======================================================================
c
c-----------------------------------------------------------------------
c  Pick an active central cell
c-----------------------------------------------------------------------
c
      if (cell_ok(ix,iy,im,jm,mask,act_val)) then
         ixwk = ix
         iywk = iy
         s=x
         t=y
        elseif (cell_ok(ix-1,iy,im,jm,mask,act_val)) then
         ixwk = ix-1
         iywk = iy
         s=c1+x
         t=y
        elseif (cell_ok(ix,iy-1,im,jm,mask,act_val)) then
         ixwk = ix
         iywk = iy-1
         s=x
         t=c1+y
        elseif (cell_ok(ix-1,iy-1,im,jm,mask,act_val)) then
         ixwk = ix-1
         iywk = iy-1
         s=c1+x
         t=c1+y
        else
         bess2d_msk = flag_val
         return
      end if
c
c-----------------------------------------------------------------------
c  Pass field and mask to 16-pt grid
c-----------------------------------------------------------------------
c
      do 20 j=1,4
         do 10 i=1,4
            fs(i,j)=f(ixwk+i-2,iywk+j-2)
            ms(i,j)=mask(ixwk+i-2,iywk+j-2)
 10      continue
 20   continue
c
c-----------------------------------------------------------------------
c  Inner points in 4x4 grid are assumed to be valid.
c  Outer points can be under the mask and are replaced by extrapolation
c  from the inner points.
c-----------------------------------------------------------------------
c
c -- Sides
c
      do 30 k=1,8
         if(ms(i0(1,k),j0(1,k)).ne.act_val) then
            fs(i0(1,k),j0(1,k))=c2*fs(i0(2,k),j0(2,k))+
     *                         cm1*fs(i0(3,k),j0(3,k))
         endif
 30   continue
c
c -- Corners
c
      do 40 k=1,4
         if(ms(i1(1,k),j1(1,k)).ne.act_val) then
            fs(i1(1,k),j1(1,k))=p5*(fs(i1(2,k),j1(2,k))+
     *                              fs(i1(3,k),j1(3,k)))
         endif
 40   continue
c
c-----------------------------------------------------------------------
c Bessel interpolation
c-----------------------------------------------------------------------
c
      bess2d_msk=bess2d(2,2,s,t,fs,4,4)
c
      return
      end