anchor.diff:      y = -x
     ELSE
       DO nd=1,nbdirs
  -      yd(nd) = 2*x*xd(nd)/(2.0*delta_x)
  -    END DO
  -    y = x**2/(2.0*delta_x) + delta_x/2.0
  +      yd(nd) = 2.0_dp*x*xd(nd)/(2.0_dp*delta_x)
  +    END DO
  +    y = x**2/(2.0_dp*delta_x) + delta_x/2.0_dp
     END IF
  -END SUBROUTINE ABS_SMOOTH_DV+END SUBROUTINE ABS_SMOOTH_DV
  +
  
positive.code:       t2d(nd) = td(nd, i+1)
        p0d(nd, i) = yptd(nd, i)
        p1d(nd, i) = t1d(nd)
        p2d(nd, i) = ((3.0*md(nd, i)-2.0*t1d(nd)-t2d(nd))*dx-(3.0*m(i)-2.0&
  &       *t1-t2)*dxd(nd))/dx**2
        p3d(nd, i) = ((t1d(nd)+t2d(nd)-2.0*md(nd, i))*dx**2-(t1+t2-2.0*m(i&
  &       ))*2*dx*dxd(nd))/(dx**2)**2
      END DO
      p0(i) = ypt(i)
      p1(i) = t1
      p2(i) = (3.0*m(i)-2.0*t1-t2)/dx
      p3(i) = (t1+t2-2.0*m(i))/dx**2
    END DO
  END SUBROUTINE SETUP_DV
  
  
  !  Differentiation of abs_smooth in forward (tangent) mode:
  !   variations   of useful results: y
  !   with respect to varying inputs: x
  SUBROUTINE ABS_SMOOTH_DV(x, xd, delta_x, y, yd, nbdirs)
  !   USE DIFFSIZES
  !  Hint: nbdirsmax should be the maximum number of differentiation directions
    IMPLICIT NONE
    INTEGER, PARAMETER :: reki=SELECTED_REAL_KIND(15, 307)
  ! in
    REAL(reki), INTENT(IN) :: x, delta_x
    REAL(reki), DIMENSION(nbdirs), INTENT(IN) :: xd
    INTEGER, INTENT(IN) :: nbdirs
  ! out
    REAL(reki), INTENT(OUT) :: y
    REAL(reki), DIMENSION(nbdirs), INTENT(OUT) :: yd
    INTRINSIC SELECTED_REAL_KIND
    INTEGER :: nd
    IF (x .GE. delta_x) THEN
      DO nd=1,nbdirs
  
negative.code:     ! evaluate the Akima spline and its derivatives
  
      implicit none
      integer, parameter :: dp = kind(0.d0)
  
      ! in
      integer, intent(in) :: npt
      real(dp), dimension(npt), intent(in) :: xpt  ! given x points
      real(dp), dimension(npt-1), intent(in) :: p0, p1, p2, p3  ! spline coefficients
      integer, intent(in) :: n
      real(dp), dimension(n), intent(in) :: x  ! x values to evalute at
      real(dp), dimension(npt-1, npt), intent(in) :: dp0dxpt, dp1dxpt, dp2dxpt, dp3dxpt
      real(dp), dimension(npt-1, npt), intent(in) :: dp0dypt, dp1dypt, dp2dypt, dp3dypt
  
      ! out
      real(dp), dimension(n), intent(out) :: y  ! interpolate y values
      real(dp), dimension(n), intent(out) :: dydx  ! derivative of y w.r.t. x
      real(dp), dimension(n, npt), intent(out) :: dydxpt, dydypt  ! derivative of y w.r.t. xpt and ypt
  
      ! local
      integer :: i, j, k
      real(dp) :: dx
  
       ! interpolate at each point
      do i = 1, n
  