! compile as: gfortran -O2 -fopenmp 
! try also the option -march=native 
 program numpi
 use omp_lib
 implicit none        
! rkind=4  single precision (32-bit)
! rkind=8  double precision (64-bit)
! rkind=10 double precision (80-bit)

 integer, parameter :: rkind=10
 integer(8), parameter:: two=2, p0=1, p1=35

 real(rkind),parameter:: zero=real(0.,rkind), half=real(0.5,rkind)
 real(rkind),parameter:: one=real(1.,rkind), four=real(4.,rkind)
 real(rkind) :: s,pi(0:1),foo,y,x,w
 integer:: ntmax,nproc,prl
 integer(8) :: m,N,power
 real(8) :: start(0:1),finish(0:1),time(0:1),tick
 character(100) :: of
 character(2):: str
! Preparations
   foo(y)=four/(one+y*y)

! Get file name for output from shell script. if not use default
   call getenv('TEST_PREC_OF',of)
!  print*, of
   if(len_trim(adjustl(of))==0) of='local_pi.dat'

! Start calculations
     write(str,'(i2)') 8*rkind
     open(1,file=of)
!       write(*,*) "# 26-digit PI= 3.1415926535897932384626433"
!       write(*,*) "# "//str//"-bit PI=acos(-1)=", acos(-one)
!       write(*,11)'# p','N=2**p  ','bit','PI (seq)','time (seq)','PI (par)','time (par)'
 
       write(1,*) "# 26-digit PI= 3.1415926535897932384626433"
!       write(1,*) "# "//str//"-bit PI=acos(-1)=", acos(-one)
       write(1,11)'# p','N=2**p  ','bit','PI (seq)','time (seq)','PI (par)','time (par)'

         do power= p0,p1
            N= two**power
            w= one/real(N,rkind)
! prl=0/1 for sequential/parallel OpenMP code
              do prl=0,1
                 s= zero
                 start(prl)=omp_get_wtime()
!!!!!!$omp parallel if(prl==1) num_threads(ntmax) shared(N,w) private(x,m) 
!$omp parallel if(prl==1) shared(N,w) private(x,m) 
!$omp do reduction(+:s)
                   do m= 1,N
                      x= w*(real(m,rkind) - half)
                      s= s + foo(x)
                   end do
!$omp end parallel
                 pi(prl)= w*s
                 finish(prl)=omp_get_wtime()
                 time(prl)=finish(prl)-start(prl)
              end do !prl
            write(1,12) power,N,rkind*8, pi(0),time(0),pi(1),time(1)
!            write(*,12) power,N,rkind*8, pi(0),time(0),pi(1),time(1)
         enddo
     close(1)

11  format(a3,a15,a4,2(a26,a15))
12  format(i3,i15,i4,2(f26.21,e15.6))
     end program