c#######################################################################
c#
c#                Mark points exceeding error tolerance.
c#   Use sum of squares of second DIFFERENCES as error indicator.
c#
c#######################################################################
      function d2error(F,Mask,errlev,ext,int,fil)
      implicit none
c-----------------------------------------------------------------------
      integer ext(4),int(4),fil(4)
      integer d2error
      real*8 F(ext(1):ext(3),ext(2):ext(4))
      real*8 Mask(ext(1):ext(3),ext(2):ext(4))
      real*8 errlev
c-----------------------------------------------------------------------
      integer i,j,imin,imax,jmin,jmax
      real*8 d2x,d2y
      real*8 ZERO,TWO
      parameter(ZERO = 0.D0,TWO = 2.D0)
c
c-----------------------------------------------------------------------
      imin = int(1)
      jmin = int(2)
      imax = int(3)
      jmax = int(4)
      if(imin.eq.fil(1)) imin = imin + 1
      if(jmin.eq.fil(2)) jmin = jmin + 1
      if(imax.eq.fil(3)) imax = imax - 1
      if(jmax.eq.fil(4)) jmax = jmax - 1
c
c-----------------------------------------------------------------------
c If patch is too small, mark all points for refinement.
      if(((fil(3)-fil(1)).lt.2).or.((fil(4)-fil(2)).lt.2)) then
        do j=int(2),int(4)
          do i=int(1),int(3)
            Mask(i,j) = ZERO
          end do
        end do
        d2error = (int(4) - int(2) + 1)*(int(3) - int(1) + 1)
        return
      end if
c
c-----------------------------------------------------------------------
      d2error = 0
      if(jmin.ne.int(2)) then
        j = jmin-1
        if(imin.ne.int(1)) then
          i = imin-1
          d2x = F(i+2,j) - TWO*F(i+1,j) + F(i,j)
          d2y = F(i,j+2) - TWO*F(i,j+1) + F(i,j)
          if((d2x*d2x + d2y*d2y).gt.errlev) then
            Mask(i,j) = ZERO
            d2error = d2error + 1
          end if
        end if
        do i = imin,imax
          d2x = F(i+1,j) - TWO*F(i,j) + F(i-1,j)
          d2y = F(i,j+2) - TWO*F(i,j+1) + F(i,j)
          if((d2x*d2x + d2y*d2y).gt.errlev) then
            Mask(i,j) = ZERO
            d2error = d2error + 1
          end if
        end do
        if(imax.ne.int(3)) then
          i = imax+1
          d2x = F(i,j) - TWO*F(i-1,j) + F(i-2,j)
          d2y = F(i,j+2) - TWO*F(i,j+1) + F(i,j)
          if((d2x*d2x + d2y*d2y).gt.errlev) then
            Mask(i,j) = ZERO
            d2error = d2error + 1
          end if
        end if
      end if
      do j = jmin,jmax
        if(imin.ne.int(1)) then
          i = imin-1
          d2x = F(i+2,j) - TWO*F(i+1,j) + F(i,j)
          d2y = F(i,j+1) - TWO*F(i,j) + F(i,j-1)
          if((d2x*d2x + d2y*d2y).gt.errlev) then
            Mask(i,j) = ZERO
            d2error = d2error + 1
          end if
        end if
        do i = imin,imax
          d2x = F(i+1,j) - TWO*F(i,j) + F(i-1,j)
          d2y = F(i,j+1) - TWO*F(i,j) + F(i,j-1)
          if((d2x*d2x + d2y*d2y).gt.errlev) then
            Mask(i,j) = ZERO
            d2error = d2error + 1
          end if
        end do
        if(imax.ne.int(3)) then
          i = imax+1
          d2x = F(i,j) - TWO*F(i-1,j) + F(i-2,j)
          d2y = F(i,j+1) - TWO*F(i,j) + F(i,j-1)
          if((d2x*d2x + d2y*d2y).gt.errlev) then
            Mask(i,j) = ZERO
            d2error = d2error + 1
          end if
        end if
      end do
      if(jmax.ne.int(4)) then
        j = jmax+1
        if(imin.ne.int(1)) then
          i = imin-1
          d2x = F(i+2,j) - TWO*F(i+1,j) + F(i,j)
          d2y = F(i,j) - TWO*F(i,j-1) + F(i,j-2)
          if((d2x*d2x + d2y*d2y).gt.errlev) then
            Mask(i,j) = ZERO
            d2error = d2error + 1
          end if
        end if
        do i = imin,imax
          d2x = F(i+1,j) - TWO*F(i,j) + F(i-1,j)
          d2y = F(i,j) - TWO*F(i,j-1) + F(i,j-2)
          if((d2x*d2x + d2y*d2y).gt.errlev) then
            Mask(i,j) = ZERO
            d2error = d2error + 1
          end if
        end do
        if(imax.ne.int(3)) then
          i = imax+1
          d2x = F(i,j) - TWO*F(i-1,j) + F(i-2,j)
          d2y = F(i,j) - TWO*F(i,j-1) + F(i,j-2)
          if((d2x*d2x + d2y*d2y).gt.errlev) then
            Mask(i,j) = ZERO
            d2error = d2error + 1
          end if
        end if
      end if
c
      return
      end
