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:           tmax=tmax*xmax
         end do
         tmpInt=zero
         if (ldebug) write(lout,"(/a,i4,6(1X,1pe16.9))") "UTILS> DEBUG polintegrate jj-step: ", &
              jj, tmax, tmin, fmax, fmin, xmax, xmin
         do kk=1,kMax-kMin+1
            tmin=tmin*xmin
            tmax=tmax*xmax
            fmin=fmin+(tmin*cof(kk))/real(kk+order-1,fPrec)
            fmax=fmax+(tmax*cof(kk))/real(kk+order-1,fPrec)
            if (ldebug) write(lout,"(a,i4,5(1X,1pe16.9))") "UTILS> DEBUG polintegrate kk-step: ", &
                 kk, tmax, tmin, fmax, fmin, cof(kk)
         end do
         select case(order)
         case(1)
  
negative.code:        end select
         polintegrate=polintegrate+tmpInt
         cumul(jj)=polintegrate
  #ifdef DEBUG
         write(lout,"(a,i4,6(1X,1pe16.9))") "UTILS> DEBUG polintegrate jj-step: ", &
              jj, tmax, tmin, fmax, fmin, tmpInt, cumul(jj)
  #endif
      end do
  
      ! fill outer range
      cumul(jMax+1:datalen)=cumul(jMax)
  
  #ifdef DEBUG
      write(lout,"(a,1pe16.9)") "UTILS> DEBUG polintegrate end: ",polintegrate
  #endif
      
    end function polintegrate
    
  