anchor.diff:  
       trace_ray: block
   
  -      real(dp), dimension(nn) :: t, lat_ray, lng_ray, h_ray, hterr, Q
  +      real(dp), dimension(nn) :: t, lat_ray, lng_ray, h_ray, hterr
         real(dp), dimension(3) :: x1, x2, xi, xw, xr
         integer :: i, j
   
  
positive.code:       call geo2xyz(src % lat, src % lng, x1(1), x1(2), x1(3))
        call geo2xyz(      lat,       lng, x2(1), x2(2), x2(3))
  
        call cross(x1, x2, xi)
        call cross(xi / sqrt(sum(xi**2)), x1, xw)
  
        interpolate_terrain: do i = 1, nn
          xi(:) = x1 * cos(t(i) * adist) + xw * sin(t(i) * adist)
  
negative.code: 
    contains
  
      pure function cross(x,y) result(z)
        real(dp), dimension(3), intent(in) :: x,y
        real(dp), dimension(3) :: z
        z(1) =   x(2)*y(3) - x(3)*y(2)
        z(2) = - x(1)*y(3) + x(3)*y(1)
        z(3) =   x(1)*y(2) - x(2)*y(1)
      end function
  
      elemental real(dp) function f(t, k)
        real(dp), intent(in) :: t, k
  