! compile as: gfortran -O2 -fopenmp summ.f90 
! try also the option -march=native 
        program summ
        use omp_lib
        implicit none        
! rkind=4  single precision (32-bit)
! rkind=8  double precision (64-bit)
! rkind=10 double precision (80-bit)
! set paral=0/1 for sequential/parallel OpenMP code
        integer, parameter :: rkind=10
        integer, parameter :: paral=0
        real(kind=rkind) :: S
        integer:: ntmax,nproc
        integer(8) :: m,N,base,p0,p1,power
        real(8) :: t0,t1,tim,tick,x
        character(100) :: of
! 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.)
! Get file name for output from shell script. if not define use default
       call getenv('SUMM_OF',of)
!       print*, of
       if(len_trim(adjustl(of))==0) of='localhost_test.dat'

! Start calculations
        base=2
        p0=10
        p1=34
        
        open(1,file=of)
        write(1,11)'# p','N=2**p','S=Sum(j,-N,N)','kind(S)','Time (s)'
!        write(*,11)'# p','N=2**p','S=Sum(j,-N,N)','kind(S)','Time (s)' !write to screen

        do power= p0,p1
           N= base**power
           S= real(0.,rkind)
           t0=omp_get_wtime()
! we do not want parallel loop if paral=0
! we use maximal possible number of parallel threads if paral=1 
!$omp parallel if(paral==1) num_threads(ntmax) private(m) shared(N) 
!$omp do reduction(+:S)
                do m= -N,N
                   S=S + real(m,rkind)
                end do
!$omp end parallel
           t1=omp_get_wtime()
           tim=t1-t0
           x=N
           write(1,12) power,x,S,kind(S),tim
!           write(*,12) power,x,S,kind(S),tim !write to screen
        enddo
        close(1)

11      format(a3,a20,a26,a8,a9)
12      format(i3,e20.12,e26.16,i4,e15.6)
        end program