anchor.diff:                                              ((2*mnuc*AtomicNumber_oper(i))/(w**2*(2*J+1)))*p_tot(w,vesc,i)
                   end if
               end do
  +            !!$OMP end do
  +            !!$OMP end parallel
           end if
       end function OMEGA_oper
   end module opermod
  
positive.code:     double precision :: ier,alist,blist,rlist,elist,iord,last!for integrator
      double precision, allocatable :: u_int_res(:)
  
      ! ! stuff for openmp testing
      ! integer :: threadnums, threadid
      ! integer, external :: omp_get_num_threads, omp_get_thread_num
  
      dimension alist(1000),blist(1000),elist(1000),iord(1000),rlist(1000)!for integrator
      !external gausstest !this is just for testing
      external integrand_oper
  
negative.code:     end do
      !$OMP end do
      !$OMP end parallel
      ! completes integral (2.4) of paper 1501.03729
      capped = 4.d0*pi*Rsun**3*capped
  
      if (capped .gt. 1.d100) then
  