function bess2d_pmsk (ix,iy,x,y,f,msk,im,jm,act_val,ncseg,mxcst,
     &                      icst,jcst,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     ICST         coast x-indices.             (integer array)      ===
c     IM, JM       1st & 2nd dimensions of F    (integer)            ===
c     JCST         coast y-indices.             (integer array)      ===
c     MSK          mask array.                  (integer array)      ===
c     MXCST        maximum length of coasts.    (integer)            ===
c     NCSEG        number of coastal segments.  (integer)            ===
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_PMSK  interpolated value           (real)               ===
c                                                                    ===
c  Calls:  BESS2D                                                    ===
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 act_val,i,im,ix,ixwk,iy,iywk,j,jm,mxcst,ncseg
      integer icst(mxcst,ncseg),jcst(mxcst,ncseg),msk(im,jm)
      FLOAT
     *      flag_val,x,y
      FLOAT
     *      fs(4,4)
      FLOAT
     *      bess2d,bess2d_pmsk
      FLOAT
     *      f(im,jm)
c
c=======================================================================
c  Begin executable code.
c=======================================================================
c
c-----------------------------------------------------------------------
c  Ensure correct constant value under extended mask.
c-----------------------------------------------------------------------
c
c     ---------------------------------------------------
c     --- Central cell is properly contained in mask. ---
c     ---------------------------------------------------
c
      iywk = min( (iy+1), jm )
      ixwk = min( (ix+1), im )
      if ( (msk(ix  ,iy  ).ne.act_val) .and.
     &     (msk(ixwk,iy  ).ne.act_val) .and.
     &     (msk(ix  ,iywk).ne.act_val) .and.
     &     (msk(ixwk,iywk).ne.act_val) ) then
         bess2d_pmsk = f(icst(1,msk(ix,iy)),jcst(1,msk(ix,iy)))
         return
      end if
c
c     ----------------------------
c     --- Point is on a coast. ---
c     ----------------------------
c
      ixwk = min( (ix+1), im )
      if ( (msk(ix  ,iy  ).ne.act_val) .and.
     &     (msk(ixwk,iy  ).ne.act_val) .and.
     &     (y.eq.c0) ) then
         bess2d_pmsk = f(icst(1,msk(ix,iy)),jcst(1,msk(ix,iy)))
         return
      end if
c
      iywk = min( (iy+1), jm )
      if ( (msk(ix  ,iy  ).ne.act_val) .and.
     &     (msk(ix  ,iywk).ne.act_val) .and.
     &     (x.eq.c0) ) then
         bess2d_pmsk = f(icst(1,msk(ix,iy)),jcst(1,msk(ix,iy)))
         return
      end if
c
c-----------------------------------------------------------------------
c  Pass field and mask to 16-pt grid
c-----------------------------------------------------------------------
c
      do 20 j=1,4
         iywk = iy+j-2
         do 10 i=1,4
            ixwk = ix+i-2
            if (msk(ixwk,iywk).eq.act_val) then
               fs(i,j)=f(ixwk,iywk)
              else
               fs(i,j) =f(icst(1,msk(ixwk,iywk)),jcst(1,msk(ixwk,iywk)))
            end if
 10      continue
 20   continue
c
c-----------------------------------------------------------------------
c Bessel interpolation
c-----------------------------------------------------------------------
c
      bess2d_pmsk = bess2d (2,2,x,y,fs,4,4)
c
      return
      end