program caf_test
  use omp_lib

  IMPLICIT NONE
  integer :: nsim
  INTEGER, PARAMETER :: K=29  ! size of array to be populated
  INTEGER :: nthreads
  INTEGER, DIMENSION(K) :: image_ind
  INTEGER :: mrem, mchunk
  INTEGER :: tid, offset
  INTEGER :: i_start, i_end
  INTEGER, ALLOCATABLE, DIMENSION(:,:) :: seq
  INTEGER :: ii, incr, n, jj, isim
  REAL(kind=8), DIMENSION(K), CODIMENSION[*] :: fill_in  ! only shared object
  INTEGER, ALLOCATABLE, DIMENSION(:) :: seed
  real(kind=8) :: tstart, tfinish
  logical :: display=.true.

  nsim = 10000

  if (display) then
    write(*,*) "Hello from image ", this_image(), "out of ", num_images()," total images"
  end if
  sync all



  nthreads = num_images()
  tid      = this_image()

  if (display.and.tid.eq.1) print *, ''

  call random_seed(size=n)
  IF (ALLOCATED(seed)) DEALLOCATE(seed)
  ALLOCATE(seed(n))
  seed = 0

  seed = 19840317*this_image() ! each image gets own seed
  call random_seed(put=seed)


  IF (ALLOCATED(seq)) DEALLOCATE(seq)
  ALLOCATE(seq(nthreads,2))

  mrem   = mod(K,nthreads)
  mchunk = (K-mrem)/nthreads

  ! create local record of image-specific indices in 'seq'
  create_seq: do ii=1,nthreads
    IF (ii-1.LT.mrem) THEN
      incr      = mchunk + 1
      offset    = 1
      seq(ii,1) = (ii-1)*incr + offset
      seq(ii,2) = seq(ii,1) + incr -1
    ELSE  ! last image does less work
      incr      = mchunk
      offset    = mrem + 1
      seq(ii,1) = (ii-1)*incr + offset
      seq(ii,2) = seq(ii,1) + incr -1
    END IF 

    if (tid.eq.1.and.display) write (*,*) 'Image ', ii, 'handles range', seq(ii,:)
    end do create_seq

  ! indices specific to each image
  if (tid-1.LT.mrem) then
    incr   = mchunk + 1
    offset = 1
  else
    incr   = mchunk
    offset = mrem + 1
  end if
  i_start  = (tid-1)*incr + offset
  i_end    = i_start + incr - 1

  sync all

  tstart = omp_get_wtime()
  busy_work: do isim=1,nsim
    fill_in = 0.0  ! 0.0 on all images
    call random_number(fill_in(i_start:i_end)) ! non-overlapping areas across images where array is populated
  
    ! synchronize coarray across images
    sync all  ! necessary
    talk: do ii=1,nthreads
      fill_in(seq(ii,1):seq(ii,2)) = fill_in(seq(ii,1):seq(ii,2))[ii] ! this is the only data transfer statement
    end do talk
    sync all

  end do busy_work

  tfinish = omp_get_wtime()

  if (display.and.tid.eq.1) then
    write (*,*), ''
    write (*,*), '----------------------------------'
    write (*,*), 'Coarray synchronized across images in', tfinish-tstart, 'seconds'
    write (*,*), '----------------------------------'
  end if

end program caf_test

