!----------------------------------------------------------------!
! The current code has the purpose:                              !
! (1) use threadprivate variable in the different program units. !
!----------------------------------------------------------------!
! Xiaoxu Guan at LSU HPC.                                        !
! March 17, 2016                                                 !
!----------------------------------------------------------------!
      module thpr
      implicit none
      integer             :: id,istart,iend
      integer, parameter  :: idp = 8, nsize = 410 
      real(kind=idp)      :: array_a(0:nsize-1),array_b(0:nsize-1),array_c(0:nsize-1)
!$omp threadprivate (id,istart,iend)
      end module thpr


      program threadprivate_v2
      use thpr
      implicit none
      integer         :: i,nothread,nchunk,omp_get_num_threads,omp_get_thread_num
      real(kind=idp)  :: diff
      real(kind=idp), dimension(0:nsize-1) :: array_master

         array_c = 0.0_idp
      do i = 1, nsize
         array_a(i) = dble(i)
         array_b(i) = dble(i+i)
      end do


!$omp parallel 
         nothread = omp_get_num_threads() 
         id = omp_get_thread_num()
         nchunk = dble(nsize)/dble(nothread)
         istart = id*nchunk
         iend = istart + nchunk - 1
      if(id == nothread-1) then
         istart = id*nchunk
         iend = nsize-1
         nchunk = iend - istart + 1
         write(*,'(1x,a,3(i3,3x))') 'Hello from ', id, istart, iend
      end if
!$omp end parallel


      write(*,'(1x,a,i3)') 'Hello I am master. id = ', id


!$omp parallel
      call scalar
!$omp end parallel


! check the results:
      do i = 0, nsize-1
         array_master(i) = array_a(i) + array_b(i)
         diff = ( array_master(i) - array_c(i) ) / array_master(i)
         if( ABS(diff) >= 1.d-10 ) then 
             write(*,'(1x,a)') 'Error seems too large!'
             write(*,'(1x,a,i4,a,es12.5,3x,2(es12.5,3x))')   &
     &            'i = ',i, ' diff = ', diff, array_master(i), array_c(i)
         end if
      end do


      end program threadprivate_v2



!--------------------------------------------------!
! this routine needs to be called by all threads.  !
!--------------------------------------------------!
      subroutine scalar
      use thpr
      implicit none
!omp threadprivate(id,istart,iend)
      integer :: i

! each thread only works on its own chunk of arrays.
      write(*,'(1x,a,3x,3(a,i3,3x))')  'from scalar rotune.',  & 
     &      "id = ", id, "istart = ", istart, "iend = ", iend

      do i = istart, iend
         array_c(i) = array_a(i) + array_b(i)
      end do

      return
      end subroutine
