anchor.diff:  
        Xout(:,p) = X_n(:)
     enddo
  +  !$acc end data
   
   end subroutine react
  
positive.code:                  else
                      I = 0.0d0
                   endif
                   
                   A(m,n) = I - dt*J(m,n)
                enddo
             enddo
             
             ! construct the RHS
             call rhs(nspec, ZERO, X_np1, dXdt, rpar, ipar)
             b(:) = X_n(:) - X_np1(:) + dt*dXdt(:)
             
             ! solve the linear system -- we'll just use Lapack.  Note that the
             ! input b is overwritten on output
             call dgesv(nspec, 1, A, nspec, ipiv, b, nspec, info)
             
             dX(:) = b(:)
             
             X_np1(:) = X_np1(:) + dX(:)
  
             ! compare to see if we converged
  
negative.code:   integer :: max_iter = 20
    integer :: iter, npts, p
    logical :: converged
  
    npts = size(T)
    ierr = 0
  
    !$acc data copyin(Xin, T, rho, tmax, eps, tol, max_iter, ic12, io16)        &
    !$acc copyout(Xout)                                                         &
    !$acc copy(ierr)
  
    !$acc kernels
  
    ierr = 0
    !$acc loop gang vector reduction(max:ierr)                                  &
    !$acc private(m,n,iter,I,time,dt,converged)                                 &
    !$acc private(dXdt, X1, X2, dX1dt, dX2dt, X_n, X_np1, dX, A, J, b)          &
    !$acc private(info, ipiv, rpar, ipar)
    do p = 1, npts
       ! get an estimate of the timestep by looking at the RHS
       ! dt = min{X/(dX/dt)}
  