!### macro's #####################################################
!
#define TRACEBACK write (gol,'("in ",a," (",a,", line",i5,")")') rname, __FILE__, __LINE__; call goErr
#define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if
#define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if
!
#include "tm5.inc"
!
!#################################################################

module OMP_ParTools

  use GO        , only : gol, goPr, goErr


  public  :: TM5_OMP_Init
  public  :: my_omp_max, my_omp_sum, my_omp_domain  
  
  ! --- const --------------------------------------

  character(len=*), parameter  ::  mname = 'OMP_ParTools'
  
  
  
  ! --- var --------------------------------------

  
  
    ! OMP variables :
   integer            :: omp_mytask, omp_ntasks, omp_mode
   logical            :: omp_master
!$OMP THREADPRIVATE (omp_mytask, omp_ntasks, omp_mode, omp_master)

   integer, parameter :: mx_omp_wrk = 32, mx_omp_ntasks = 64
   real               :: wrk(mx_omp_wrk,0:1,0:mx_omp_ntasks-1)

  
  
contains


  ! ===================================================
  
  ! -----------------------------------------------------------------  
! OMP subroutines, originally written by r. richter, sgi, july 2008
! adapted for TM5-structure and conventions by Vincent Huijnen, september 2008
! -----------------------------------------------------------------

   subroutine TM5_OMP_Init ( status )

! Written : r. richter, sgi, july 2008

   implicit none

!vh #if defined (MPI)
!vh #include <mpif.h>
!vh #endif

    
    ! --- in/out ----------------------------------
    
    integer, intent(out)            ::  status
    
    ! --- const ------------------------------
    
! Local parameter

   integer            :: i, j, mype, npes, info

! Functions

   integer            :: omp_get_thread_num, omp_get_num_threads
   logical            :: omp_in_parallel

!vh #if defined (_OPENMP)
#ifdef _OPENMP
   if (omp_in_parallel()) &
      stop 'CALL my_omp_init () in parrallel region is not allowed'
#endif

!vh#if defined (MPI)
!vh   call mpi_comm_rank (mpi_comm_world, mype, info)
!vh   call mpi_comm_size (mpi_comm_world, npes, info)

!vh   call mpi_barrier (mpi_comm_world, info)
!vh   do i = 0, npes - 1
!vh      if (i == mype) then
!vh#endif

!$OMP PARALLEL &
!$OMP  default (none) &
!$OMP  shared  (wrk) &
!$OMP  private (j)

#ifdef _OPENMP
         j = omp_get_thread_num()
         omp_mytask = omp_get_thread_num()
         omp_ntasks = omp_get_num_threads()
         write(*,*)'OMP. ntasks: ',omp_ntasks,', my_task: ',omp_mytask
#else
         j = 0
         omp_mytask = 0
         omp_ntasks = 1
         write(*,*)'No openmp parallelization used'
#endif

         if (omp_ntasks > mx_omp_ntasks) then
            print *, 'ERROR : Increase mx_omp_ntasks to', omp_ntasks
            stop
         end if

         omp_mode   = 0
         omp_master = omp_mytask == 0

         wrk(1:mx_omp_wrk,0:1,j) = 0.0

!$OMP END PARALLEL

!vh#if defined (MPI)
!vh      end if
!vh      call mpi_barrier (mpi_comm_world, info)
!vh   end do
!vh#endif
    ! ok
    status = 0
   return
   end subroutine TM5_OMP_Init

! 

   subroutine my_omp_domain (std_start, std_end, dom_start, dom_end)

! Written : r. richter, sgi, july 2008

   implicit none

! Input/output parameters

   integer            :: std_start, std_end, dom_start, dom_end

! Local parameters

   real (8)           :: x_len, x_mytask, x_ntasks

! Functions

   logical            :: omp_in_parallel


!vh #if defined (_OPENMP)
#ifdef _OPENMP
   if (omp_in_parallel()) then
      x_mytask = real (omp_mytask, 8)
      x_ntasks = real (omp_ntasks, 8)
   else
      x_mytask = 0.0_8
      x_ntasks = 1.0_8
   end if
#else
   x_mytask = real (omp_mytask, 8)
   x_ntasks = real (omp_ntasks, 8)
#endif

   x_len    = real (std_end - std_start + 1, 8)

   dom_start = std_start + floor ((x_mytask * x_len + 0.5_8) / &
               x_ntasks, 4)
   dom_end   = std_start - 1 + floor ((x_mytask * x_len + x_len + &
               0.5_8) / x_ntasks, 4)

   return
   end subroutine my_omp_domain

! 

   subroutine my_omp_sum (lvec, gvec)

! Written : r. richter, sgi, july 2008

   implicit none

! Input/output parameters

   real               :: lvec, gvec

! Local parameters

   logical            :: do_parallel

! Functions

   logical (4)        :: omp_in_parallel


!#if defined (_OPENMP)
#ifdef _OPENMP
   do_parallel = omp_in_parallel()
#else
   do_parallel = .false.
#endif

   if (do_parallel) then

      omp_mode = mod (omp_mode + 1, 2)
      wrk(1,omp_mode,omp_mytask) = lvec

!$OMP BARRIER

      if (omp_master) gvec = sum (wrk(1,omp_mode,0:omp_ntasks-1))

!$OMP BARRIER

   else
      gvec = lvec
   end if

   return
   end subroutine my_omp_sum

! 

   subroutine my_omp_max (lvec, gvec)

! Written : r. richter, sgi, july 2008

   implicit none

! Input/output parameters

   real               :: lvec, gvec

! Local parameters

   logical            :: do_parallel

! Functions

   logical (4)        :: omp_in_parallel


!#if defined (_OPENMP)
#ifdef _OPENMP
   do_parallel = omp_in_parallel()
#else
   do_parallel = .false.
#endif

   if (do_parallel) then

      omp_mode = mod (omp_mode + 1, 2)
      wrk(1,omp_mode,omp_mytask) = lvec

!$OMP BARRIER

      if (omp_master) gvec = maxval (wrk(1,omp_mode,0:omp_ntasks-1))

!$OMP BARRIER

   else
      gvec = lvec
   end if

   return
   end subroutine my_omp_max


end module OMP_ParTools