function spl2d2(xb,yb,nbrx,x,nbry,y,f,fx,fy,fxy,maxy,ixd,iyd)
c-----------------------------------------------------------------------
      dimension p(5,4),f(maxy,nbrx),fx(maxy,nbrx),fy(maxy,nbrx),
     1 fxy(maxy,nbrx),z(3),t(4),u(4),x(nbrx),y(nbry)
      integer i1,ix,iy,i2,j1,j2
      real c,d,al2,bl2,xl2,bra,brb
c-----------------------------------------------------------------------
      data p/3*0.,-1.,4*0.,1.,4*0.,-1.,4*0.,1.,0./
      z(1)=xb
      z(3)=yb
      if (z(1).le.x(1)) go to 1
      if (z(1).ge. x(nbrx)) go to 2
      call search (z(1),x,nbrx,j,miss)
   50 if (z(3).le.y(1)) go to 3
      if (z(3).ge.y(nbry)) go to 4
      call search (z(3),y,nbry,i,miss)
    5 c=1./6.
      u(1)=x(j)
      u(3)=y(i)
      d=1./3.
      u(2)=x(j+1)-u(1)
      u(4)=y(i+1)-u(3)
      do 10 i1=1,3,2
      p(3,i1+1)=z(i1)-u(i1)
      p(3,i1)=u(i1+1)-p(3,i1+1)
      al2=p(3,i1+1)**2
      bl2=p(3,i1)**2
      xl2=u(i1+1)**2
      bra=(al2-xl2)*c
      brb=(bl2-xl2)*c
      p(1,i1+1)=bra*p(3,i1+1)
      p(1,i1)=brb*p(3,i1)
      p(2,i1+1)=bra+al2*d
   10 p(2,i1)=-brb-bl2*d
      ix=ixd+1
      iy=iyd+1
      do 15 i1=1,4
   15 t(i1)=0.
      do 20 i1=1,2
      i2=i+i1-1
      do 20 j1=1,2
      j2=j1+j-1
      t(i1)=t(i1)+fxy(i2,j2)*p(ix,j1)
      t(i1)=t(i1)+fy(i2,j2)*p(ix+2,j1)
      t(i1+2)=t(i1+2)+fx(i2,j2)*p(ix,j1)
   20 t(i1+2)=t(i1+2)+f(i2,j2)*p(ix+2,j1)
      spl2d2=(t(1)*p(iy,3)+t(2)*p(iy,4)+t(3)*p(iy+2,3)
     1       +t(4)*p(iy+2,4))/(u(2)*u(4))
      return
c---take care of extrapolation
    1 j=1
      go to 50
    2 j=nbrx-1
      go to 50
    3 i=1
      go to 5
    4 i=nbry-1
      go to 5
      end
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++