anchor.diff:        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
  +      p2d(nd, i) = ((3.0_dp*md(nd, i)-2.0_dp*t1d(nd)-t2d(nd))*dx-(3.0_dp&
  +&       *m(i)-2.0_dp*t1-t2)*dxd(nd))/dx**2
  +      p3d(nd, i) = ((t1d(nd)+t2d(nd)-2.0_dp*md(nd, i))*dx**2-(t1+t2-&
  +&       2.0_dp*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
  +    p2(i) = (3.0_dp*m(i)-2.0_dp*t1-t2)/dx
  +    p3(i) = (t1+t2-2.0_dp*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)
  +  INTEGER, PARAMETER :: dp=KIND(0.d0)
   ! in
  -  REAL(reki), INTENT(IN) :: x, delta_x
  -  REAL(reki), DIMENSION(nbdirs), INTENT(IN) :: xd
  -  INTEGER, INTENT(IN) :: nbdirs
  +  REAL(dp), INTENT(IN) :: x, delta_x
  +  REAL(dp), 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
  +  REAL(dp), INTENT(OUT) :: y
  +  REAL(dp), DIMENSION(nbdirs), INTENT(OUT) :: yd
  +  INTRINSIC KIND
     INTEGER :: nd
     IF (x .GE. delta_x) THEN
       DO nd=1,nbdirs
  
positive.code:     elseif (x <= -delta_x) then
          y = -x
      else
          y = x**2/(2.0*delta_x) + delta_x/2.0
      end if
  
  end subroutine abs_smooth
  
negative.code: 
      implicit none
  
      integer, parameter :: dp = kind(0.d0)
      real(dp), parameter :: eps = 1d-30
  
      ! in
      integer, intent(in) :: n
      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(dp), dimension(n-1), intent(out) :: p0, p1, p2, p3  ! spline coefficients
  
      ! local
      integer :: i
      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
  