anchor.diff:    deallocate (output_array)
     deallocate (res_reference)
   
  -!  reference values for self test
  +  ! reference values for self test
     allocate (res_reference(num_variables))
     res_reference(1) = -2.0795476275887927d0
     res_reference(2) = 2.1489169824673575d-2
  
positive.code:     res = output_array(vector_length, 1)
      print *, 'dE/dx_i for i =', ideriv, 'is', res
  
  !     test against reference numbers
      if (abs(res - res_reference(ideriv)) > 1.0d-6) then
        print *, 'error: derivatives do not match reference numbers'
        stop 1
  
negative.code:   endif
  
    vector_length = 2**order
    allocate (density(vector_length, num_variables, num_points))
    density = 0.0d0
    do ipoint = 1, num_points
      ! we use fantasy values here
      density(1, 1, ipoint) = 1.0d0 !        n zeroth order
      density(1, 2, ipoint) = 2.0d0 !nabla_x n zeroth order
      density(1, 3, ipoint) = 3.0d0 !nabla_y n zeroth order
  