! 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

 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,foo,y,x,w

 integer:: ntmax,nproc,ip
 integer(8), parameter:: itwo=2
 integer(8) :: m,N,power,p0,p1
 real(8) :: start,tick
 real:: time
 character(100) :: of
 character(2):: str

 !!! pi26=3.1415926535897932384626433_16
! Preparations
   foo(y)=four/(one+y*y)

! Set OpenMP options
!   ntmax=omp_get_max_threads()
!   nproc=omp_get_num_procs()
!   tick=omp_get_wtick()
!   call omp_set_dynamic(.false.)
!   call omp_set_nested(.false.)

! set the file name for output
   write(str,'(i2)') 8*rkind
   of="PI_"//str//".dat"
! check if the file exists. if it does not, create it.   
   call system( "f="//trim(adjustl(of))//"; [ -f $f ] || touch $f" )

!!!!!!! Start calculations
     open(1,file=of,position='append') ! do not overwrite file but rather append new records
!     open(1,file=of) !overwrite file

         p0=31
         p1=40

         do power= p0,p1
            N= itwo**power
            w= one/real(N,rkind)
            s= zero
            start=omp_get_wtime()
!!!!!!!!!!!$omp parallel num_threads(ntmax) shared(N,w) private(x,m) 
!$omp parallel 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= w*s
            time=omp_get_wtime()-start
            ip=power
            write(1,*) ip,N,pi !,time 
            write(*,*) ip,N,pi !,time 
         enddo
     close(1)
     
     end program