! compile as: gfortran -O2 -fopenmp sumj-omp.f 
! try also the option -march=native 
        program sumjomp
        use omp_lib
        implicit none        
        integer, parameter :: Prec=10, prl=0
! Prec=8  double precision (64-bit)
! Prec=10 double precision (80-bit)
! prl=0 sequential code
! prl=1 use OpenMP for parallel calculation
        integer:: ntmax,nproc
        integer(8) ::  j,N,base,p0,p1,power
        real(8) :: t0,t1,t,tick
        real(Prec) ::  x
        character(3) :: sp
        character(100) :: of,h
! 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 output file 
       write(sp,'(I3)') Prec
       call HOSTNM(h)
       if(prl==0) then
         of=trim(adjustl(sp))//"_"//trim(adjustl(h))//".dat"
       elseif(prl==1) then
         of=trim(adjustl(sp))//"_"//trim(adjustl(h))//"_omp.dat"
       else
         print*,"Illegal value of prl=",prl
         stop
       end if
! Collect system info and write that to file
       call system("bash sys_info.sh "//of)

!!!!! Actual calculation
        
        base=2
        p0=10
        p1=34
        open(1,file=trim(of),position="append")
        do power= p0,p1
           N= base**power
           x= real(0.,Prec)
           t0=omp_get_wtime()
! we do not want parallel loop if prl=0
!$omp parallel num_threads(ntmax) if(prl==1) private(j) shared(N) 
!$omp do reduction(+:x)
                do j= -N,N
                   x=x + real(j,Prec)
                end do
!$omp end parallel
           t1=omp_get_wtime()
           t=t1-t0
           write(1,*) power,N,x,kind(x),t
!           write(*,*) power,N,x,kind(x),t
        enddo 
        close(1)
        end program