anchor.diff:    integer :: i
   
     ! Formatting rules
  +  character(*), parameter :: format1='(i3,1(a,e13.6))'
     character(*), parameter :: format3='(i3,3(a,e13.6))'
     character(*), parameter :: format4='(i3,4(a,e13.6))'
     character(*), parameter :: format6='(i3,6(a,e13.6))'
  
positive.code:   real(RP) :: dataout9(9, 10)  ! results from Velocity Gradient
    real(RP) :: dataout18(18, 10) ! results from Velocity Hessian
  
  
    ! Declare the return type of the turblib functions as integer.
    ! This is required for custom error handling (see the README).
    integer :: getvelocity, getforce, getvelocityandpressure, getvelocitygradient
    integer :: getboxfilter, getboxfiltersgs, getboxfiltergradient
    integer :: getvelocitylaplacian, getvelocityhessian
    integer :: getpressuregradient, getpressurehessian
    integer :: getposition
    ! If working with hdf5 cutout files, uncomment the lines below
    ! to load the file and use the functions locally:
  
negative.code:        ', d2pdxdz=', dataout6(3,i), ', d2pdydy=', dataout6(4,i),  &
         ', d2pdydz=', dataout6(5,i), ', d2pdzdz', dataout6(6,i)
    end do
  
   
    ! getrawvelocity
    !write(*,*)
    !write(*,'(a)') 'Requesting raw velocity at 10 points...'
    !rc = getrawvelocity(authkey, dataset,  time, x, y, z, xwidth, ywidth, zwidth, rawvelocity)
    !do i = 1, 10
    !  write(*,format3) i, ': Vx=', rawvelocity(3*i), ', Vy=', rawvelocity(3*i+1), ', Vz=', rawvelocity(3*i+2)
    !end do
  
    ! getrawpressure
  
    write(*,*)
    write(*,'(2(a,f8.6),a)') 'Requesting position at 10 points, starting at time ', &
  