! compile as: gfortran -O2 -fopenmp 
! try also the option -march=native 
 program calce
 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=8

 real(rkind),parameter:: one=real(1.,rkind) 
 real(rkind) :: x,w
 
 integer :: ntmax,nproc,ip
 real :: time
 integer(8), parameter:: itwo=2
 integer(8) :: m,N,power,p0,p1
 real(8) :: start,finish,tick
 character(100) :: of
 character(2):: str
 
 !!! e26=2.7182818284590452353602874_16
 
! Preparations

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

     end program