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(6),int(6),fil(6)
      integer d2error
      real*8 F(ext(1):ext(4),ext(2):ext(5),ext(3):ext(6))
      real*8 Mask(ext(1):ext(4),ext(2):ext(5),ext(3):ext(6))
      real*8 errlev
c-----------------------------------------------------------------------
      integer k,i,j,kmin,kmax,imin,imax,jmin,jmax
      real*8 d2w,d2x,d2y
      real*8 ZERO,TWO
      parameter(ZERO = 0.D0,TWO = 2.D0)
c
c-----------------------------------------------------------------------
      kmin = int(1)
      imin = int(2)
      jmin = int(3)
      kmax = int(4)
      imax = int(5)
      jmax = int(6)
      if(kmin.eq.fil(1)) kmin = kmin + 1
      if(imin.eq.fil(2)) imin = imin + 1
      if(jmin.eq.fil(3)) jmin = jmin + 1
      if(kmax.eq.fil(4)) kmax = kmax - 1
      if(imax.eq.fil(5)) imax = imax - 1
      if(jmax.eq.fil(6)) jmax = jmax - 1
c
c-----------------------------------------------------------------------
c If patch is too small, mark all points for refinement.
      if(((fil(4)-fil(1)).lt.2).or.((fil(5)-fil(2)).lt.2).or.
     .                             ((fil(6)-fil(3)).lt.2)) then
        do j=int(3),int(6)
          do i=int(2),int(5)
            do k=int(1),int(4)
              Mask(k,i,j) = ZERO
            end do
          end do
        end do
        d2error = (int(6) - int(3) + 1)*(int(5) - int(2) + 1)
     .                                 *(int(4) - int(1) + 1)
        return
      end if
c
c-----------------------------------------------------------------------
      d2error = 0
      if(jmin.ne.int(3)) then
        j = jmin-1
        if(imin.ne.int(2)) then
          i = imin-1
          if(kmin.ne.int(1)) then
            k = kmin-1
            d2w = F(k+2,i,j) - TWO*F(k+1,i,j) + F(k,i,j)
            d2x = F(k,i+2,j) - TWO*F(k,i+1,j) + F(k,i,j)
            d2y = F(k,i,j+2) - TWO*F(k,i,j+1) + F(k,i,j)
            if((d2w*d2w + d2x*d2x + d2y*d2y).gt.errlev) then
              Mask(k,i,j) = ZERO
              d2error = d2error + 1
            end if
          end if
          do k = kmin,kmax
            d2w = F(k+1,i,j) - TWO*F(k,i,j) + F(k-1,i,j)
            d2x = F(k,i+2,j) - TWO*F(k,i+1,j) + F(k,i,j)
            d2y = F(k,i,j+2) - TWO*F(k,i,j+1) + F(k,i,j)
            if((d2w*d2w + d2x*d2x + d2y*d2y).gt.errlev) then
              Mask(k,i,j) = ZERO
              d2error = d2error + 1
            end if
          end do
          if(kmax.ne.int(4)) then
            k = kmax+1
            d2w = F(k,i,j) - TWO*F(k-1,i,j) + F(k-2,i,j)
            d2x = F(k,i+2,j) - TWO*F(k,i+1,j) + F(k,i,j)
            d2y = F(k,i,j+2) - TWO*F(k,i,j+1) + F(k,i,j)
            if((d2w*d2w + d2x*d2x + d2y*d2y).gt.errlev) then
              Mask(k,i,j) = ZERO
              d2error = d2error + 1
            end if
          end if
        end if
        do i = imin,imax
          if(kmin.ne.int(1)) then
            k = kmin-1
            d2w = F(k+2,i,j) - TWO*F(k+1,i,j) + F(k,i,j)
            d2x = F(k,i+1,j) - TWO*F(k,i,j) + F(k,i-1,j)
            d2y = F(k,i,j+2) - TWO*F(k,i,j+1) + F(k,i,j)
            if((d2w*d2w + d2x*d2x + d2y*d2y).gt.errlev) then
              Mask(k,i,j) = ZERO
              d2error = d2error + 1
            end if
          end if
          do k = kmin,kmax
            d2w = F(k+1,i,j) - TWO*F(k,i,j) + F(k-1,i,j)
            d2x = F(k,i+1,j) - TWO*F(k,i,j) + F(k,i-1,j)
            d2y = F(k,i,j+2) - TWO*F(k,i,j+1) + F(k,i,j)
            if((d2w*d2w + d2x*d2x + d2y*d2y).gt.errlev) then
              Mask(k,i,j) = ZERO
              d2error = d2error + 1
            end if
          end do
          if(kmax.ne.int(4)) then
            k = kmax+1
            d2w = F(k,i,j) - TWO*F(k-1,i,j) + F(k-2,i,j)
            d2x = F(k,i+1,j) - TWO*F(k,i,j) + F(k,i-1,j)
            d2y = F(k,i,j+2) - TWO*F(k,i,j+1) + F(k,i,j)
            if((d2w*d2w + d2x*d2x + d2y*d2y).gt.errlev) then
              Mask(k,i,j) = ZERO
              d2error = d2error + 1
            end if
          end if
        end do
        if(imax.ne.int(5)) then
          i = imax+1
          if(kmin.ne.int(1)) then
            k = kmin-1
            d2w = F(k+2,i,j) - TWO*F(k+1,i,j) + F(k,i,j)
            d2x = F(k,i,j) - TWO*F(k,i-1,j) + F(k,i-2,j)
            d2y = F(k,i,j+2) - TWO*F(k,i,j+1) + F(k,i,j)
            if((d2w*d2w + d2x*d2x + d2y*d2y).gt.errlev) then
              Mask(k,i,j) = ZERO
              d2error = d2error + 1
            end if
          end if
          do k = kmin,kmax
            d2w = F(k+1,i,j) - TWO*F(k,i,j) + F(k-1,i,j)
            d2x = F(k,i,j) - TWO*F(k,i-1,j) + F(k,i-2,j)
            d2y = F(k,i,j+2) - TWO*F(k,i,j+1) + F(k,i,j)
            if((d2w*d2w + d2x*d2x + d2y*d2y).gt.errlev) then
              Mask(k,i,j) = ZERO
              d2error = d2error + 1
            end if
          end do
          if(kmax.ne.int(4)) then
            k = kmax+1
            d2w = F(k,i,j) - TWO*F(k-1,i,j) + F(k-2,i,j)
            d2x = F(k,i,j) - TWO*F(k,i-1,j) + F(k,i-2,j)
            d2y = F(k,i,j+2) - TWO*F(k,i,j+1) + F(k,i,j)
            if((d2w*d2w + d2x*d2x + d2y*d2y).gt.errlev) then
              Mask(k,i,j) = ZERO
              d2error = d2error + 1
            end if
          end if
        end if
      end if
      do j = jmin,jmax
        if(imin.ne.int(2)) then
          i = imin-1
          if(kmin.ne.int(1)) then
            k = kmin-1
            d2w = F(k+2,i,j) - TWO*F(k+1,i,j) + F(k,i,j)
            d2x = F(k,i+2,j) - TWO*F(k,i+1,j) + F(k,i,j)
            d2y = F(k,i,j+1) - TWO*F(k,i,j) + F(k,i,j-1)
            if((d2w*d2w + d2x*d2x + d2y*d2y).gt.errlev) then
              Mask(k,i,j) = ZERO
              d2error = d2error + 1
            end if
          end if
          do k = kmin,kmax
            d2w = F(k+1,i,j) - TWO*F(k,i,j) + F(k-1,i,j)
            d2x = F(k,i+2,j) - TWO*F(k,i+1,j) + F(k,i,j)
            d2y = F(k,i,j+1) - TWO*F(k,i,j) + F(k,i,j-1)
            if((d2w*d2w + d2x*d2x + d2y*d2y).gt.errlev) then
              Mask(k,i,j) = ZERO
              d2error = d2error + 1
            end if
          end do
          if(kmax.ne.int(4)) then
            k = kmax+1
            d2w = F(k,i,j) - TWO*F(k-1,i,j) + F(k-2,i,j)
            d2x = F(k,i+2,j) - TWO*F(k,i+1,j) + F(k,i,j)
            d2y = F(k,i,j+1) - TWO*F(k,i,j) + F(k,i,j-1)
            if((d2w*d2w + d2x*d2x + d2y*d2y).gt.errlev) then
              Mask(k,i,j) = ZERO
              d2error = d2error + 1
            end if
          end if
        end if
        do i = imin,imax
          if(kmin.ne.int(1)) then
            k = kmin-1
            d2w = F(k+2,i,j) - TWO*F(k+1,i,j) + F(k,i,j)
            d2x = F(k,i+1,j) - TWO*F(k,i,j) + F(k,i-1,j)
            d2y = F(k,i,j+1) - TWO*F(k,i,j) + F(k,i,j-1)
            if((d2w*d2w + d2x*d2x + d2y*d2y).gt.errlev) then
              Mask(k,i,j) = ZERO
              d2error = d2error + 1
            end if
          end if
          do k = kmin,kmax
            d2w = F(k+1,i,j) - TWO*F(k,i,j) + F(k-1,i,j)
            d2x = F(k,i+1,j) - TWO*F(k,i,j) + F(k,i-1,j)
            d2y = F(k,i,j+1) - TWO*F(k,i,j) + F(k,i,j-1)
            if((d2w*d2w + d2x*d2x + d2y*d2y).gt.errlev) then
              Mask(k,i,j) = ZERO
              d2error = d2error + 1
            end if
          end do
          if(kmax.ne.int(4)) then
            k = kmax+1
            d2w = F(k,i,j) - TWO*F(k-1,i,j) + F(k-2,i,j)
            d2x = F(k,i+1,j) - TWO*F(k,i,j) + F(k,i-1,j)
            d2y = F(k,i,j+1) - TWO*F(k,i,j) + F(k,i,j-1)
            if((d2w*d2w + d2x*d2x + d2y*d2y).gt.errlev) then
              Mask(k,i,j) = ZERO
              d2error = d2error + 1
            end if
          end if
        end do
        if(imax.ne.int(5)) then
          i = imax+1
          if(kmin.ne.int(1)) then
            k = kmin-1
            d2w = F(k+2,i,j) - TWO*F(k+1,i,j) + F(k,i,j)
            d2x = F(k,i,j) - TWO*F(k,i-1,j) + F(k,i-2,j)
            d2y = F(k,i,j+1) - TWO*F(k,i,j) + F(k,i,j-1)
            if((d2w*d2w + d2x*d2x + d2y*d2y).gt.errlev) then
              Mask(k,i,j) = ZERO
              d2error = d2error + 1
            end if
          end if
          do k = kmin,kmax
            d2w = F(k+1,i,j) - TWO*F(k,i,j) + F(k-1,i,j)
            d2x = F(k,i,j) - TWO*F(k,i-1,j) + F(k,i-2,j)
            d2y = F(k,i,j+1) - TWO*F(k,i,j) + F(k,i,j-1)
            if((d2w*d2w + d2x*d2x + d2y*d2y).gt.errlev) then
              Mask(k,i,j) = ZERO
              d2error = d2error + 1
            end if
          end do
          if(kmax.ne.int(4)) then
            k = kmax+1
            d2w = F(k,i,j) - TWO*F(k-1,i,j) + F(k-2,i,j)
            d2x = F(k,i,j) - TWO*F(k,i-1,j) + F(k,i-2,j)
            d2y = F(k,i,j+1) - TWO*F(k,i,j) + F(k,i,j-1)
            if((d2w*d2w + d2x*d2x + d2y*d2y).gt.errlev) then
              Mask(k,i,j) = ZERO
              d2error = d2error + 1
            end if
          end if
        end if
      end do
      if(jmax.ne.int(6)) then
        j = jmax+1
        if(imin.ne.int(2)) then
          i = imin-1
          if(kmin.ne.int(1)) then
            k = kmin-1
            d2w = F(k+2,i,j) - TWO*F(k+1,i,j) + F(k,i,j)
            d2x = F(k,i+2,j) - TWO*F(k,i+1,j) + F(k,i,j)
            d2y = F(k,i,j) - TWO*F(k,i,j-1) + F(k,i,j-2)
            if((d2w*d2w + d2x*d2x + d2y*d2y).gt.errlev) then
              Mask(k,i,j) = ZERO
              d2error = d2error + 1
            end if
          end if
          do k = kmin,kmax
            d2w = F(k+1,i,j) - TWO*F(k,i,j) + F(k-1,i,j)
            d2x = F(k,i+2,j) - TWO*F(k,i+1,j) + F(k,i,j)
            d2y = F(k,i,j) - TWO*F(k,i,j-1) + F(k,i,j-2)
            if((d2w*d2w + d2x*d2x + d2y*d2y).gt.errlev) then
              Mask(k,i,j) = ZERO
              d2error = d2error + 1
            end if
          end do
          if(kmax.ne.int(4)) then
            k = kmax+1
            d2w = F(k,i,j) - TWO*F(k-1,i,j) + F(k-2,i,j)
            d2x = F(k,i+2,j) - TWO*F(k,i+1,j) + F(k,i,j)
            d2y = F(k,i,j) - TWO*F(k,i,j-1) + F(k,i,j-2)
            if((d2w*d2w + d2x*d2x + d2y*d2y).gt.errlev) then
              Mask(k,i,j) = ZERO
              d2error = d2error + 1
            end if
          end if
        end if
        do i = imin,imax
          if(kmin.ne.int(1)) then
            k = kmin-1
            d2w = F(k+2,i,j) - TWO*F(k+1,i,j) + F(k,i,j)
            d2x = F(k,i+1,j) - TWO*F(k,i,j) + F(k,i-1,j)
            d2y = F(k,i,j) - TWO*F(k,i,j-1) + F(k,i,j-2)
            if((d2w*d2w + d2x*d2x + d2y*d2y).gt.errlev) then
              Mask(k,i,j) = ZERO
              d2error = d2error + 1
            end if
          end if
          do k = kmin,kmax
            d2w = F(k+1,i,j) - TWO*F(k,i,j) + F(k-1,i,j)
            d2x = F(k,i+1,j) - TWO*F(k,i,j) + F(k,i-1,j)
            d2y = F(k,i,j) - TWO*F(k,i,j-1) + F(k,i,j-2)
            if((d2w*d2w + d2x*d2x + d2y*d2y).gt.errlev) then
              Mask(k,i,j) = ZERO
              d2error = d2error + 1
            end if
          end do
          if(kmax.ne.int(4)) then
            k = kmax+1
            d2w = F(k,i,j) - TWO*F(k-1,i,j) + F(k-2,i,j)
            d2x = F(k,i+1,j) - TWO*F(k,i,j) + F(k,i-1,j)
            d2y = F(k,i,j) - TWO*F(k,i,j-1) + F(k,i,j-2)
            if((d2w*d2w + d2x*d2x + d2y*d2y).gt.errlev) then
              Mask(k,i,j) = ZERO
              d2error = d2error + 1
            end if
          end if
        end do
        if(imax.ne.int(5)) then
          i = imax+1
          if(kmin.ne.int(1)) then
            k = kmin-1
            d2w = F(k+2,i,j) - TWO*F(k+1,i,j) + F(k,i,j)
            d2x = F(k,i,j) - TWO*F(k,i-1,j) + F(k,i-2,j)
            d2y = F(k,i,j) - TWO*F(k,i,j-1) + F(k,i,j-2)
            if((d2w*d2w + d2x*d2x + d2y*d2y).gt.errlev) then
              Mask(k,i,j) = ZERO
              d2error = d2error + 1
            end if
          end if
          do k = kmin,kmax
            d2w = F(k+1,i,j) - TWO*F(k,i,j) + F(k-1,i,j)
            d2x = F(k,i,j) - TWO*F(k,i-1,j) + F(k,i-2,j)
            d2y = F(k,i,j) - TWO*F(k,i,j-1) + F(k,i,j-2)
            if((d2w*d2w + d2x*d2x + d2y*d2y).gt.errlev) then
              Mask(k,i,j) = ZERO
              d2error = d2error + 1
            end if
          end do
          if(kmax.ne.int(4)) then
            k = kmax+1
            d2w = F(k,i,j) - TWO*F(k-1,i,j) + F(k-2,i,j)
            d2x = F(k,i,j) - TWO*F(k,i-1,j) + F(k,i-2,j)
            d2y = F(k,i,j) - TWO*F(k,i,j-1) + F(k,i,j-2)
            if((d2w*d2w + d2x*d2x + d2y*d2y).gt.errlev) then
              Mask(k,i,j) = ZERO
              d2error = d2error + 1
            end if
          end if
        end if
      end if
c
      return
      end
