anchor.diff:      ! actually interpolate
       polinterp = nevilleMethod(xvals(jMin:jMax),yvals(jMin:jMax),mpoints,x,dy)
   
  -    if ( ldebug ) then
  -      write(lout,'(/a,1pe16.9)')   'UTILS> interpolating arrays for x= ',x
  -      write(lout,'(a,i0)')         'UTILS>    arrays of length: ',datalen
  -      write(lout,'(a,i0)')         'UTILS>    initial guess: jguess=',jguess
  -      write(lout,'(a,i0)')         'UTILS>    index found at:  jj=',jj
  -      write(lout,'(2(a,1pe16.9))') 'UTILS>    i.e. between x= ',xvals(jj),' and x= ',xvals(jj+1)
  -      write(lout,'(3(a,i0))')      'UTILS>    using ',mpoints,' points: jMin= ',jMin,' and jMax= ',jMax
  -      write(lout,'(2(a,1pe16.9))') 'UTILS>    i.e. between x= ',xvals(jMin),'and x= ',xvals(jMax)
  -      write(lout,'(a,1pe16.9)')    'UTILS> calculated value: ',polinterp
  -      write(lout,'(a)')            'UTILS> arrays:'
  -      do ii=1,datalen
  -        write(lout,'(i0,2(1X,1pe16.9))') ii, xvals(ii), yvals(ii)
  -      enddo
  -    end if
  +#ifdef DEBUG
  +    write(lout,'(/a,1pe16.9)')   'UTILS> interpolating arrays for x= ',x
  +    write(lout,'(a,i0)')         'UTILS>    arrays of length: ',datalen
  +    write(lout,'(a,i0)')         'UTILS>    initial guess: jguess=',jguess
  +    write(lout,'(a,i0)')         'UTILS>    index found at:  jj=',jj
  +    write(lout,'(2(a,1pe16.9))') 'UTILS>    i.e. between x= ',xvals(jj),' and x= ',xvals(jj+1)
  +    write(lout,'(3(a,i0))')      'UTILS>    using ',mpoints,' points: jMin= ',jMin,' and jMax= ',jMax
  +    write(lout,'(2(a,1pe16.9))') 'UTILS>    i.e. between x= ',xvals(jMin),'and x= ',xvals(jMax)
  +    write(lout,'(a,1pe16.9)')    'UTILS> calculated value: ',polinterp
  +    write(lout,'(a)')            'UTILS> arrays:'
  +    do ii=1,datalen
  +      write(lout,'(i0,2(1X,1pe16.9))') ii, xvals(ii), yvals(ii)
  +    enddo
  +#endif
       
       ! update guess for next call
       jguess=jj
  
positive.code:     do ii=1, datalen-1
        if(xvals(ii) >= xvals(ii+1)) then
          checkArray=.false.
          if (llprint) write(lerr,"(a)") "UTILS> ERROR checkArray: xvals should be in increasing order"
        end if
      end do
    end function checkArray
  
negative.code:     jMax=datalen
      if(present(rmax)) jMax=huntBin(rmax,xvals,datalen,-1)
  
  #ifdef DEBUG
      write(lout,"(/a,5(1X,i0))") "UTILS> DEBUG polintegrate start: ",datalen,mpoints,order,jMin,jMax
  #endif
  
      do jj=jMin+1,jMax ! loop over data points
         kMin=min(max(jj-1-(mpoints-1)/2,1),datalen+1-mpoints)
  