anchor.diff:    real(RP) :: dataout9(9, 10)  ! results from Velocity Gradient
     real(RP) :: dataout18(18, 10) ! results from Velocity Hessian
   
  +  integer,parameter :: x=0, y=0, z=0, xwidth=16, ywidth=16, zwidth=16
  +  real(RP),dimension(xwidth*ywidth*zwidth*3) :: rawvelocity 
   
     ! 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 :: getpressure, getpressuregradient, getpressurehessian
  +  integer :: getrawvelocity, getrawpressure
     integer :: getposition
     ! If working with hdf5 cutout files, uncomment the lines below
     ! to load the file and use the functions locally:
  
positive.code:        ', d2pdxdz=', dataout6(3,i), ', d2pdydy=', dataout6(4,i),  &
         ', d2pdydz=', dataout6(5,i), ', d2pdzdz', dataout6(6,i)
    end do
  
    write(*,*)
    write(*,'(2(a,f8.6),a)') 'Requesting position at 10 points, starting at time ', &
  
negative.code:   end do
  
    write(*,*)
    write(*,'(a)') 'Requesting pressure at 10 points...'
    rc = getpressure(authkey, dataset,  time, Lag6, NoTInt, 10, points, dataout1)
    do i = 1, 10 
      write(*,format1) i, ': ', dataout1(i)
    end do
  
    write(*,*)
    write(*,'(a)') 'Requesting pressure gradient at 10 points...'
    rc = getpressuregradient(authkey, dataset,  time, FD4Lag4, NoTInt, 10, points, dataout3)
    do i = 1, 10 
  