anchor.diff:    ! the density variable of interest to 1, and set other perturbed
     ! densities to 0
   
  -  allocate(density(vector_length, num_variables, num_points))
  -  density = 0.0d0
  -  do ipoint = 1, num_points
  +  allocate(density(vector_length, num_density_variables, num_grid_points))
  +  density = 0.0d0
  +  do ipoint = 1, num_grid_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
  
positive.code:   end do
  
    res = derivative(id, &
                     num_points, &
                     num_variables, &
                     vector_length, &
                     density)
    deallocate(density)
  
negative.code:   end do
  
    res = derivative(id, &
                     num_grid_points, &
                     num_density_variables, &
                     vector_length, &
                     density)
    deallocate(density)
  