anchor.diff:      TSpline1D_Derivative = dely/ho+ ((1-3*a0**2)*this%ddF(llo) + (3*b0**2-1)*this%ddF(lhi))*ho/6
   
       end function TSpline1D_Derivative
  +
  +    subroutine TCubicSpline_IntegralArray(this, yint, first_index, last_index)
  +    !Array of values of integral at the original sampled points
  +    class(TCubicSpline) :: this
  +    real(sp_acc), intent(out) :: yint(*)
  +    integer, intent(in), optional :: first_index, last_index
  +    integer first, last
  +    real(sp_acc) dx
  +    integer i, ix
  +
  +    First  = PresentDefault(1, first_index)
  +    Last = PresentDefault(this%n, last_index)
  +    yint(1) = 0
  +    ix = 2
  +    do i=first+1, last
  +        dx = (this%x(i) - this%x(i-1))
  +        yint(ix) = yint(ix-1) + dx*( (this%F(i)+this%F(i-1))/2 - dx**2/24*(this%ddF(i)+this%ddF(i-1)))
  +        ix = ix+1
  +    end do
  +
  +    end subroutine TCubicSpline_IntegralArray
   
       subroutine TSpline1D_ArrayValue(this, x, y, error )
       !Get array of values y(x), assuming x is monotonically increasing
  
positive.code:     procedure :: FindValue => TCubicSpline_FindValue
      procedure :: LoadState => TCubicSpline_LoadState
      procedure :: SaveState => TCubicSpline_SaveState
      generic :: Init => TCubicSpline_Init, TCubicSpline_InitInt
      FINAL :: TCubicSpline_Free !not needed in standard, just for compiler bugs
      end Type
  
negative.code:             ! direction but outside in the y direction.
          ELSE IF ((ixdi > 0.AND.ixdi < nxd) .AND.  &
              (iydi <= 0.OR.iydi >= nyd)) THEN
              ! Retrieves the z and partial derivative values at the other
              ! vertex of the semi-infinite rectangle.
              IF (ixdi /= ixdipv .OR. iydi /= iydipv) THEN
                  ixd1 = ixd0 + 1
                  dx = xd(ixd1) - x0
                  dxsq = dx*dx
                  z10 = zd(ixd1,iyd0)
                  zx10 = pdd(1,ixd1,iyd0)
                  zy10 = pdd(2,ixd1,iyd0)
                  zxy10 = pdd(3,ixd1,iyd0)
                  ! Calculates the polynomial coefficients.
                  z0dx = (z10-z00)/dx
                  zy0dx = (zy10-zy00)/dx
                  p00 = z00
                  p01 = zy00
                  p10 = zx00
                  p11 = zxy00
                  p20 = (2.0* (z0dx-zx00)+z0dx-zx10)/dx
                  p21 = (2.0* (zy0dx-zxy00)+zy0dx-zxy10)/dx
                  p30 = (-2.0*z0dx+zx10+zx00)/dxsq
                  p31 = (-2.0*zy0dx+zxy10+zxy00)/dxsq
              END IF
              ! Evaluates the polynomial.
              u = xii - x0
              v = yii - y0
              q0 = p00 + v*p01
              q1 = p10 + v*p11
              q2 = p20 + v*p21
              q3 = p30 + v*p31
              zii = q0 + u* (q1+u* (q2+u*q3))
              ! End of Case 2
  
              ! Case 3.  When the rectangle is outside the data area in the x
              ! direction but inside in the y direction.
          ELSE IF ((ixdi <= 0.OR.ixdi >= nxd) .AND.  &
              (iydi > 0 .AND. iydi < nyd)) THEN
              ! Retrieves the z and partial derivative values at the other
              ! vertex of the semi-infinite rectangle.
              IF (ixdi /= ixdipv .OR. iydi /= iydipv) THEN
                  iyd1 = iyd0 + 1
                  dy = yd(iyd1) - y0
                  dysq = dy*dy
                  z01 = zd(ixd0,iyd1)
                  zx01 = pdd(1,ixd0,iyd1)
                  zy01 = pdd(2,ixd0,iyd1)
                  zxy01 = pdd(3,ixd0,iyd1)
                  ! Calculates the polynomial coefficients.
                  z0dy = (z01-z00)/dy
                  zx0dy = (zx01-zx00)/dy
                  p00 = z00
                  p01 = zy00
                  p02 = (2.0*(z0dy-zy00)+z0dy-zy01)/dy
                  p03 = (-2.0*z0dy+zy01+zy00)/dysq
                  p10 = zx00
                  p11 = zxy00
                  p12 = (2.0*(zx0dy-zxy00) + zx0dy - zxy01)/dy
                  p13 = (-2.0*zx0dy + zxy01 + zxy00)/dysq
              END IF
  
              ! Evaluates the polynomial.
              u = xii - x0
              v = yii - y0
              q0 = p00 + v* (p01 + v*(p02+v*p03))
              q1 = p10 + v* (p11 + v*(p12+v*p13))
              zii = q0 + u*q1
              ! End of Case 3
  
              ! Case 4.  When the rectangle is outside the data area in both the
              ! x and y direction.
          ELSE IF ((ixdi <= 0 .OR. ixdi >= nxd) .AND.  &
              (iydi <= 0 .OR. iydi >= nyd)) THEN
              ! Calculates the polynomial coefficients.
              IF (ixdi /= ixdipv .OR. iydi /= iydipv) THEN
                  p00 = z00
                  p01 = zy00
                  p10 = zx00
                  p11 = zxy00
              END IF
              ! Evaluates the polynomial.
              u = xii - x0
              v = yii - y0
              q0 = p00 + v*p01
              q1 = p10 + v*p11
              zii = q0 + u*q1
          END IF
          ! End of Case 4
          zi(iip) = zii
  