! 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
 integer(8), parameter:: two=2, p0=1, p1=35
! integer(8), parameter:: two=2, p0=37, p1=40

 real(rkind),parameter:: one=real(1.,rkind) 
 real(rkind) :: e(0:1),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

! 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_e.dat'

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

       write(1,*) "# 26-digit E= 2.7182818284590452353602874" 
!       write(1,*) "# "//str//"-bit exp(1)=", exp(one)
       write(1,11)'# p','N = 2**p  ','bit','E (seq)','time(seq)','E (par)','time(par)'
 
         do power= p0,p1
            N= two**power
            w= one + one/real(N,rkind)
              do prl=0,1 ! prl=0/1 for sequential/parallel OpenMP code
                 x= one
                 start(prl)=omp_get_wtime()
!!!!!$omp parallel if(prl==1) num_threads(nth) shared(N,w) private(m) 
!$omp parallel if(prl==1) shared(N,w) private(m) 
!$omp do reduction(*:x)
                   do m= 1,N
                      x= x * w
                   end do
!$omp end parallel
                 finish(prl)=omp_get_wtime()
                 time(prl)=finish(prl) - start(prl)
                 e(prl)= x
              end do !prl
!              write(*,12) power,N,rkind*8,e(0),time(0),e(1),time(1)
              write(1,12) power,N,rkind*8,e(0),time(0),e(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