subroutine bssl_interp(fld_in,undef,gxout,gyout,
     $     nii,nji,nio,njo,fld_out,iunit_diag,
     $     cyclicxi,spole_point,npole_point,bessel,istat)

      dimension fld_in(nii,nji),fld_out(nio,njo),
     $     gxout(nio+1),gyout(njo+1)

      dimension fr(4)

      logical cyclicxi,diag_out,bessel,spole_point,npole_point
      
      diag_out=.false.
ccc      diag_out=.true.
      jchk=2
      istat=1
C         
C         convert the box boundaries to grid point center
C         
      do i=1,nio
C         
C         check for crossing the boundaries 
C         the only way for this to occur is if the field is 
C         cyclically continuous in x
C

        if(gxout(i+1).lt.gxout(i)) then
          gxout(i)=((gxout(i)-float(nii))+gxout(i+1))*0.5+0.5
          if(gxout(i).lt.1.0) gxout(i)=float(nii)+gxout(i)
        else
          gxout(i)=(gxout(i)+gxout(i+1))*0.5+0.5
          if(gxout(i).lt.1.0) gxout(i)=float(nii)+gxout(i)
        endif
      end do

      do j=1,njo
        gyout(j)=(gyout(j)+gyout(j+1))*0.5+0.5
C         
C         check if a pole points on the input grid
C         
        if(spole_point.and.gyout(j).lt.1.0) gyout(j)=1.0
        if(npole_point.and.gyout(j).gt.nji+0.5) gyout(j)=nji

      end do

      ibb=1
      iee=nio
      jbb=1
      jee=njo

      niim1=nii-1
      njim1=nji-1
      
      do j=jbb,jee
        do i=ibb,iee
          
          ic=ifix(gxout(i))
          jc=ifix(gyout(j))
          icp1=ic+1
          if(cyclicxi.and.icp1.gt.nii) icp1=icp1-nii 
          jcp1=jc+1
          
          if(diag_out.and.j.eq.jchk) write(iunit_diag,*)
     $         'i,j,gxout,gyout,ic,jc',
     $         i,j,gxout(i),gyout(j),ic,jc
          
          if((jc.lt.1.or.jc.gt.nji).or.
     $         (.not.cyclicxi.and.(ic.lt.1.or.ic.gt.nii))) then
            fld_out(i,j)=undef
            if(diag_out.and.j.eq.jchk) 
     $           write(iunit_diag,*) 'out of bounds ' 
            go to 100
          end if
C         
C------------------------------------------------
C         
C         bilinear/bessel interpolation based on the FNOC routine
C         bssl5 by D. Hensen, FNOC
C         
C------------------------------------------------
C         
          r=gxout(i)-ic
          s=gyout(j)-jc
C         
C         interior check
C         
          if((jc.ge.2.and.jc.lt.njim1.and.cyclicxi).or.
     $         (ic.ge.2.and.jc.ge.2.and.
     $         ic.lt.niim1.and.jc.lt.njim1)) go to 10
C         
C         border zone check
C         
          if((jc.lt.nji.and.cyclicxi).or.
     $         (ic.lt.nii.and.jc.lt.nji)) go to 5 
C         
C------------------------------------------------
C         
C         top and right edge processing
C         
C------------------------------------------------
C         
          if(ic.eq.nii.and.jc.eq.nji) then
            
            fld_out(i,j)=fld_in(nii,nji)
            if(diag_out.and.j.eq.jchk) 
     $           write(iunit_diag,*) 'upper right hand corenr'
            
          else if(ic.eq.nii) then
            
            if(diag_out.and.j.eq.jchk) 
     $           write(iunit_diag,*) 'right edge'
            
            if(fld_in(ic,jc).ne.undef.and.
     $           fld_in(ic,jcp1).ne.undef) then
              fld_out(i,j)=(1.0-s)*fld_in(ic,jc)+
     $             s*fld_in(ic,jcp1)
            else
              fld_out(i,j)=undef
            endif
            
          else if(jc.eq.nji) then
            
            if(diag_out.and.j.eq.jchk) 
     $           write(iunit_diag,*) 'top edge'
            
            if(fld_in(ic,jc).ne.undef.and.
     $           fld_in(icp1,jc).ne.undef) then
              fld_out(i,j)=(1.0-r)*fld_in(ic,jc)+
     $             r*fld_in(icp1,jc)
            else
              fld_out(i,j)=undef
            endif
            
          endif
          
          go to 100
          
 5        continue
C         
C------------------------------------------------
C         
C         border zone; bilinear
C         
C------------------------------------------------
C         
          iok_bilinear=1
          if(fld_in(ic,jc).eq.undef.or.
     $         fld_in(icp1,jc).eq.undef.or.
     $         fld_in(ic,jcp1).eq.undef.or.
     $         fld_in(icp1,jcp1).eq.undef) iok_bilinear=0
          
          if(diag_out.and.j.eq.jchk) 
     $         write(iunit_diag,*) 
     $         'border zone, iok_bilinear= ',iok_bilinear
                    
          if(iok_bilinear.eq.0) then
            fld_out(i,j)=undef
          else
            fld_out(i,j)=(1.-s)*((1.-r)*fld_in(ic,jc)
     $           +r*fld_in(icp1,jc))
     $           +s*((1.-r)*fld_in(ic,jcp1)
     $           +r*fld_in(icp1,jcp1))
          endif

          go to 100
 
 10       continue
C         
C------------------------------------------------
C         
C         interior zone
C         
C------------------------------------------------
C         
C         first check if bilinear is OK
C         
          iok_bilinear=1
          if(fld_in(ic,jc).eq.undef.or.
     $         fld_in(icp1,jc).eq.undef.or.
     $         fld_in(ic,jcp1).eq.undef.or.
     $         fld_in(icp1,jcp1).eq.undef) iok_bilinear=0
          
          if(diag_out.and.j.eq.jchk) 
     $         write(iunit_diag,*) 
     $         'interior zone, iok_bilinear= ',iok_bilinear
          
          if(iok_bilinear.eq.0) then
            fld_out(i,j)=undef
            go to 100
          else 
C         
C         bilinear value is the first guess
C         
            fld_out(i,j)=(1.-s)*((1.-r)*fld_in(ic,jc)
     $           +r*fld_in(icp1,jc))
     $           +s*((1.-r)*fld_in(ic,jcp1)
     $           +r*fld_in(icp1,jcp1))
C         
C         exit if only doing bilinear
C         
            if(.not.bessel) go to 100
            
          endif
          
          
C         
C         interpolate 4 columns (i-1,i,i+1,i+2) to j+s and store in fr(1)
C         through fr(4)
C         
          r1=r-0.5
          r2=r*(r-1.)*0.5
          r3=r1*r2*0.3333333333334
          s1=s-0.5
          s2=s*(s-1.)*0.5
          s3=s1*s2*0.3333333333334
C         
          k=0
          im1=ic-1
          ip2=ic+2
          
          if(diag_out.and.j.eq.jck) write(iunit_diag,*) 'bessel interp'
          
C         
          do ii=im1,ip2
            
            k=k+1
            
            i1=ii
            
            if(cyclicxi.and.i1.lt.1) i1=nii-i1
            if(cyclicxi.and.i1.gt.nii) i1=i1-nii
            
            j1=jc
            
            j1p1=j1+1
            j1p2=j1+2
            j1m1=j1-1
            
            fij=fld_in(i1,j1)
            fijp1=fld_in(i1,j1p1)
            fijp2=fld_in(i1,j1p2)
            fijm1=fld_in(i1,j1m1)
            if(diag_out.and.j.eq.jchk) 
     $           write(iunit_diag,*) 'i1,j1 = ',i1,j1,
     $           fij,fijp1,fijp2,fijm1
C         
C         exit if any value undefined
C         
            if(fij.eq.undef.or.
     $           fijp1.eq.undef.or.
     $           fijp2.eq.undef.or.
     $           fijm1.eq.undef) go to 100 
            
            u=(fij+fijp1)*0.5
            del=fijp1-fij
            udel2=(fijp2-fijp1+fijm1-fij)*0.5
            del3=fijp2-fijp1-2.*del+fij-fijm1
            
            fr(k)=u+s1*del+s2*udel2+s3*del3
            
          end do
          
C         
C         interpolate the fr row to ii+r
C         
          u=(fr(2)+fr(3))*0.5
          del=fr(3)-fr(2)
          udel2=(fr(4)-fr(3)+fr(1)-fr(2))*0.5
          del3=fr(4)-fr(3)-2.*del+fr(2)-fr(1)
          
          fld_out(i,j)=u+r1*del+r2*udel2+r3*del3
          
          
 100      continue
          if(diag_out.and.j.eq.jchk) 
     $         write(iunit_diag,*) 'interp value = ',fld_out(i,j)
          
        end do
      end do
      
 999  continue
      
      return
      end