   #pragma unroll
   for (int j = 0; j < NSP; ++j) {
     double yj_orig = y[j];
-    double r = fmax(srur * fabs(yj_orig), r0 / error[j]);
+    double r = fmax(srur * fabs(yj_orig), r0 / ewt[j]);
     
     #if FD_ORD == 1
       y[j] = yj_orig + r;

  
  #pragma unroll
  for (int i = 0; i < NSP; ++i) {
    error[i] = ATOL + (RTOL * fabs(cy[i]));
    y[i] = cy[i];
  }


  double sum = 0.0;
  #pragma unroll
  for (int i = 0; i < NSP; ++i) {
    sum += (ewt[i] * dy[i]) * (ewt[i] * dy[i]);
  }
  double fac = sqrt(sum / ((double)(NSP)));
  double r0 = 1000.0 * RTOL * DBL_EPSILON * ((double)(NSP)) * fac;
