anchor.diff:    real(RP) :: dataout9(9, 10)  ! results from Gradient function
     real(RP) :: dataout18(18, 10) ! results from Hessian function
   
  +  integer,parameter :: x=0, y=0, z=0, xwidth=16, ywidth=16, zwidth=16
  +  !real(RP) :: rawdata(xwidth*ywidth*zwidth*3)
  +  real(RP) :: rawdata(3,xwidth*ywidth*zwidth)
  +  real(RP) :: rawpressure(xwidth*ywidth*zwidth)
   
     ! Declare the return type of the turblib functions as integer.
     ! This is required for custom error handling (see the README).
  
positive.code: 
    ! 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))'
    character(*), parameter :: format9='(i3,9(a,e13.6))'
  
negative.code:   end do
  
    write(*,*)
    write(*,'(a)') 'Requesting raw velocity ...'
    rc = getrawvelocity(authkey, dataset,  time, x, y, z, xwidth, ywidth, zwidth, rawdata)
    do i = 1, xwidth*ywidth*zwidth
      !write(*,rawformat3) i, ': Vx=', rawdata(3*(i-1)+1), ', Vy=', rawdata(3*(i-1)+2), ', Vz=', rawdata(3*(i-1)+3)
      !write(*,rawformat3) i, ': Vx=', rawdata(1,i), ', Vy=', rawdata(2,i), ', Vz=', rawdata(3,i)
    end do
  
    write(*,*)
    write(*,'(a)') 'Requesting raw pressure ...'
    rc = getrawpressure(authkey, dataset,  time, x, y, z, xwidth, ywidth, zwidth, rawpressure)
    do i = 1, xwidth*ywidth*zwidth
      !write(*,rawformat1) i, ': ', rawpressure(i)
    end do
  
    write(*,*)
    write(*,'(a)') 'Requesting raw magnetic field ...'
    rc = getrawmagneticfield(authkey, dataset,  time, x, y, z, xwidth, ywidth, zwidth, rawdata)
    do i = 1, xwidth*ywidth*zwidth
      !write(*,rawformat3) i, ': Vx=', rawdata(3*(i-1)+1), ', Vy=', rawdata(3*(i-1)+2), ', Vz=', rawdata(3*(i-1)+3)
      !write(*,rawformat3) i, ': Vx=', rawdata(1,i), ', Vy=', rawdata(2,i), ', Vz=', rawdata(3,i)
    end do
  
    write(*,*)
    write(*,'(a)') 'Requesting raw vector potential ...'
    rc = getrawvectorpotential(authkey, dataset,  time, x, y, z, xwidth, ywidth, zwidth, rawdata)
    do i = 1, xwidth*ywidth*zwidth
      !write(*,rawformat3) i, ': Vx=', rawdata(3*(i-1)+1), ', Vy=', rawdata(3*(i-1)+2), ', Vz=', rawdata(3*(i-1)+3)
      !write(*,rawformat3) i, ': Vx=', rawdata(1,i), ', Vy=', rawdata(2,i), ', Vz=', rawdata(3,i)
    end do
  
    write(*,*)
    write(*,'(2(a,f8.6),a)') 'Requesting position at 10 points, starting at time ', &
         startTime, ' and ending at time ', endTime, '...'
    rc = getposition(authkey, dataset, startTime, endTime, lag_dt, Lag6, 10, points, dataout3)
  