anchor.diff:  
       implicit none
   
  -    integer, parameter :: ReKi = selected_real_kind(15, 307)
  -    real(ReKi), parameter :: eps = 1d-30
  +    integer, parameter :: dp = kind(0.d0)
  +    real(dp), parameter :: eps = 1d-30
   
       ! in
       integer, intent(in) :: n
  -    real(ReKi), dimension(n), intent(in) :: xpt, ypt  ! given points
  -    real(ReKi), intent(in) :: delta_x
  -    !f2py real(ReKi), intent(in) :: delta_x = 0.1
  +    real(dp), dimension(n), intent(in) :: xpt, ypt  ! given points
  +    real(dp), intent(in) :: delta_x
  +    !f2py real(dp), intent(in) :: delta_x = 0.1
   
       ! out
  -    real(ReKi), dimension(n-1), intent(out) :: p0, p1, p2, p3  ! spline coefficients
  +    real(dp), dimension(n-1), intent(out) :: p0, p1, p2, p3  ! spline coefficients
   
       ! local
       integer :: i
  -    real(ReKi), dimension(-1:n+1) :: m
  -    real(ReKi), dimension(n) :: t
  -    real(ReKi) :: m1, m2, m3, m4, w1, w2
  -    real(ReKi) :: t1, t2, dx
  +    real(dp), dimension(-1:n+1) :: m
  +    real(dp), dimension(n) :: t
  +    real(dp) :: m1, m2, m3, m4, w1, w2
  +    real(dp) :: t1, t2, dx
   
       ! compute segment slopes
       do i = 1, n-1
  
positive.code: 
  
  !        Generated by TAPENADE     (INRIA, Tropics team)
  !  Tapenade 3.8 (r5001) -  5 Nov 2013 10:56
  !
  !  Differentiation of setup in forward (tangent) mode:
  !   variations   of useful results: p0 p1 p2 p3
  !   with respect to varying inputs: xpt ypt
  !   RW status of diff variables: xpt:in p0:out p1:out p2:out p3:out
  !                ypt:in
  SUBROUTINE SETUP_DV(n, xpt, xptd, ypt, yptd, p0, p0d, p1, p1d, p2, p2d, &
  & p3, p3d, delta_x, nbdirs)
  !   USE DIFFSIZES
  !  Hint: nbdirsmax should be the maximum number of differentiation directions
    IMPLICIT NONE
    INTEGER, PARAMETER :: reki=SELECTED_REAL_KIND(15, 307)
    REAL(reki), PARAMETER :: eps=1d-30
  ! in
    INTEGER, INTENT(IN) :: n, nbdirs
  ! given points
    REAL(reki), DIMENSION(n), INTENT(IN) :: xpt, ypt
    REAL(reki), DIMENSION(nbdirs, n), INTENT(IN) :: xptd, yptd
    REAL(reki), INTENT(IN) :: delta_x
  !f2py real(ReKi), intent(in) :: delta_x = 0.1
  ! out
  ! spline coefficients
    REAL(reki), DIMENSION(n - 1), INTENT(OUT) :: p0, p1, p2, p3
    REAL(reki), DIMENSION(nbdirs, n-1), INTENT(OUT) :: p0d, p1d, p2d, &
  & p3d
  ! local
    INTEGER :: i
    REAL(reki), DIMENSION(-1:n + 1) :: m
    REAL(reki), DIMENSION(nbdirs, -1:n+1) :: md
    REAL(reki), DIMENSION(n) :: t
    REAL(reki), DIMENSION(nbdirs, n) :: td
    REAL(reki) :: m1, m2, m3, m4, w1, w2
    REAL(reki), DIMENSION(nbdirs) :: m1d, m2d, m3d, m4d, w1d, w2d
    REAL(reki) :: t1, t2, dx
    REAL(reki), DIMENSION(nbdirs) :: t1d, t2d, dxd
    INTRINSIC SELECTED_REAL_KIND
    REAL(reki) :: arg1
    REAL(reki), DIMENSION(nbdirs) :: arg1d
    INTEGER :: nd
  !   INTEGER :: nbdirs
    DO nd=1,nbdirs
      md(nd, :) = 0.0
    END DO
  ! compute segment slopes
    DO i=1,n-1
  
negative.code:   END DO
    DO nd=1,nbdirs
  ! estimation for end points
      md(nd, 0) = 2.0_dp*md(nd, 1) - md(nd, 2)
      md(nd, -1) = 2.0_dp*md(nd, 0) - md(nd, 1)
      md(nd, n) = 2.0_dp*md(nd, n-1) - md(nd, n-2)
      md(nd, n+1) = 2.0_dp*md(nd, n) - md(nd, n-1)
    END DO
    m(0) = 2.0_dp*m(1) - m(2)
    m(-1) = 2.0_dp*m(0) - m(1)
    m(n) = 2.0_dp*m(n-1) - m(n-2)
    m(n+1) = 2.0_dp*m(n) - m(n-1)
    DO nd=1,nbdirs
      td(nd, :) = 0.0_dp
    END DO
  ! slope at points
    DO i=1,n
  