anchor.diff:    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)}
  
positive.code: 
          ! initial guess
          X_np1(:) = X_n(:)
       
          do iter = 1, max_iter
  
             ! construct a numerical Jacobian via simple differencing
             ! Lapack matrix ordering is described here:
             ! http://www.netlib.org/lapack/lug/node116.html#subsecarrayargs
             ! in general, A(i,j) holds matrix element a_{i,j}
          
             ! J_{m,n} = df_m/dX_n
             do n = 1, nspec
                X1(:) = X_np1(:)
                X1(n) = X1(n) + eps
                
                X2(:) = X_np1(:)
                X2(n) = X2(n) - eps
                
                call rhs(nspec, ZERO, X1, dX1dt, rpar, ipar)
                call rhs(nspec, ZERO, X2, dX2dt, rpar, ipar)
  
  
negative.code: 
       dt = 1.d33
       dt = 0.01*min(dt, abs(Xin(ih1,p)/dXdt(ih1)))
  
       ! do the integration
       time = ZERO
       do while (time < tmax)
  