anchor.diff:  
     ! Formatting rules
     character(*), parameter :: format1='(i3,1(a,e13.6))'
  +  character(*), parameter :: rawformat1='(i4,1(a,e13.6))'
     character(*), parameter :: format3='(i3,3(a,e13.6))'
  +  character(*), parameter :: rawformat3='(i4,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))'
  
positive.code:   integer :: getvectorpotential, getvectorpotentialgradient
    integer :: getvectorpotentiallaplacian, getvectorpotentialhessian
    integer :: getpressuregradient, getpressurehessian
    integer :: getposition
    integer :: getpressure
    ! If working with hdf5 cutout files, uncomment the lines below
  
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)
  