anchor.diff:  
        Xout(:,p) = X_n(:)
     enddo
  +  !$acc end data
   
   end subroutine react
  
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)
  