# Modern Fortran Quadpack

**URL:** https://fortran-lang.discourse.group/t/modern-fortran-quadpack/3714
**Category:** Announcements
**Created:** [June 10, 2022, 3:22pm UTC](https://fortran-lang.discourse.group/t/modern-fortran-quadpack/3714 "2022-06-10T15:22:32Z")
**Posts on this page:** 1
**Showing post:** 16

<div class="post-metadata">

### Author: ![sshestov](https://avatars.discourse-cdn.com/v4/letter/s/bcef8e/32.png) [@sshestov](https://fortran-lang.discourse.group/u/sshestov)
#### Post date: [September 17, 2026, 2:22pm UTC](https://fortran-lang.discourse.group/t/modern-fortran-quadpack/3714/16 "2026-09-17T14:22:41Z")

</div>

Dear all, I’m having troubles with OpenMP-parallelized do-loops with `QAWO` or `QAGS`. My minimal (not working) example is below

```fortran

  !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(r, r_Aplane, int_cos, int_sin, f_cos, f_sin, val, eps, neval, ierr, last, iwork, work)
  do j=1, 500  
    r=A_axis(j) !! A-plane coord.
    !! ************ my straightforward sampling of the 1st plane works well with OpenMP
    f_cos = 2. * dpi * r_axis * cos(dpi*r_axis**2 *lamz0_1) * BESSEL_J0(2.*dpi*r*r_axis*lamz0_1)
    f_sin = 2. * dpi * r_axis * sin(dpi*r_axis**2 *lamz0_1) * BESSEL_J0(2.*dpi*r*r_axis*lamz0_1)
    int_cos = integrate_trap(real(r_axis,kind=8), f_cos)      
    int_sin = integrate_trap(real(r_axis,kind=8), f_sin)

    !! ********but I want to implement it via QUAD_PACK routines**************************
    !! **************variable substitution in the integral ksi <=> r_axis*****************
    !! Int ( ksi*cos( .. ksi^2) * J0(2pi r ksi/lamz0) dksi ) =>
    !! 1/2 Int ( cos(ksi^2) * J0(2pi r sqrt(ksi^2) / lamz0) dksi^2
    r_Aplane=r
    call dqags(qags_fc,0._wp,710.0_wp,1d-2, 1d-4 ,val, eps, neval, ierr, limit, lenw, last, iwork, work)
    int_cos = val*0.5*2. ! 2 to conform the definition in my original code 
    
    r_Aplane=r
    call dqags(qags_fs,0._wp,710.0_wp,1d-2, 1d-4 ,val, eps, neval, ierr, limit, lenw, last, iwork, work)
    int_sin = val*0.5*2.
    
    !!!! temporarily put newly computed real part to imaginary or vise versa
    !!psi_imag(j) = 1. - ( cos(dpi*r**2*lamz0_1)*int_sin + sin(dpi*r**2*lamz0_1)*int_cos ) * lamz0_1  
    psi_real(j) = ( cos(dpi*r**2*lamz0_1)*int_cos - sin(dpi*r**2*lamz0_1)*int_sin ) * lamz0_1
    if (mod(j,100) == 0) write(*,'(I4,A)', advance='no') int(j/100)," "
  end do
  !$OMP END PARALLEL DO
...
  contains
  real(wp) function qags_fc(x)
    real(wp), intent(in) :: x
    qags_fc = 2. * dpi * x * cos(dpi * x**2 * lamz0_1) * BESSEL_J0(2.*dpi*r_Aplane*x*lamz0_1)
  end function qags_fc

end

```

When I use self-implemented trapezoid implementation and explicit sampling (in r\_axis) it works flawlessly. But I wanted to switch to QUAD\_PACK routines or a separate Levin method, with a hope that presence of cos/sin terms can speed up.

I’ve started from Levin method - it was working but giving spurious outliers; then QAWO was working only until some j. I guess this is because Bessel J0 is also oscillating. Finally I switched to QAGS, which seems to work, but I struggling with OpenMP.

It works, but gives completely wrong results.

Do you have any ideas/suggestions?

PS: this is to calculate diffraction behind a circular occulter.

---

_[View the full topic](https://fortran-lang.discourse.group/t/modern-fortran-quadpack/3714)._