!-------------------------------------------------------------------- 
!                    Test  parallelization of random numbers 
!-------------------------------------------------------------------- 
Module Rands
Contains
!-------------------------------------------               random number generator
      FUNCTION RANDd(M)
!------------------------------------------------
      DATA LC,AM,KI,K1,K2,K3,K4,L1,L2,L3,L4  &
       /453815927,2147483648.,2147483647,536870912,131072,256, &
        16777216,4,16384,8388608,128/
      ML=M/K1*K1
      M1=(M-ML)*L1
      ML=M/K2*K2
      M2=(M-ML)*L2
      ML=M/K3*K3
      M3=(M-ML)*L3
      ML=M/K4*K4
      M4=(M-ML)*L4
      M5=KI-M
      IF(M1.GE.M5)M1=M1-KI-1
      ML=M+M1
      M5=KI-ML
      IF(M2.GE.M5)M2=M2-KI-1
      ML=ML+M2
      M5=KI-ML
      IF(M3.GE.M5)M3=M3-KI-1
      ML=ML+M3
      M5=KI-ML
      IF(M4.GE.M5)M4=M4-KI-1
      ML=ML+M4
      M5=KI-ML
      IF(LC.GE.M5)ML=ML-KI-1
      M=ML+LC
      RANDd=M/AM
      RETURN
      END FUNCTION RANDd
end Module Rands
!-------------------------------------------------------------------------------
Module Halos
      PARAMETER (Nproc = 16 ) ! total number of cores
      PARAMETER (Nhalos = 320 ) ! total number of halos
      PARAMETER (Nselect = 140 ) ! number of selected particles
      PARAMETER (Nskip    = 1e5 ) ! number of seeds to skip

      DIMENSION NumSeed(Nproc),NhaloSeeds(Nselect,Nhalos)
      Real :: Fraction
Contains
SUBROUTINE Halo(iHalo,iproc)
   use Rands
      iseed  = NumSeed(iproc)
      write(*,'(i5,3x,i3,3x,i11," i=",i4," Rand=",g12.4)') iHalo,iproc,iseed
   Do i =1, Nselect
      xx =Randd(iseed)
      If(xx<Fraction)NhaloSeeds(i,iHalo) =0
      write(*,'(i5,3x,i3,3x,i11," i=",i4," Rand=",g12.4)') iHalo,iproc,iseed,i,xx
   EndDo
   NumSeed(iproc) = iseed
end SUBROUTINE Halo
end Module Halos
!
!-------------------------------------------------------------------------------------
Program RandTest
  use Rands
  use Halos

  NhaloSeeds =1
      iseed = 9121071
      Do i=1,Nproc
         NumSeed(i) = iseed
         Do j=1,Nskip
            xx = Randd(iseed)
         EndDo
      EndDo
      write(*,'(4i12)') NumSeed
      Fraction =0.10
      LengthChunk =INT(Nhalos/Nproc)
!$OMP PARALLEL DO DEFAULT(SHARED) &    
!$OMP PRIVATE (ic,iproc)
Do ic =1,Nhalos 
   iproc  =INT((ic-1)/LengthChunk)+1
   CALL  Halo(ic,iproc)
EndDo

Do i=1,Nhalos
   isum =SUM(NhaloSeeds(:,i))
   write(*,'(i5,f8.4,3x,200i1)') i,(1.-float(isum)/Nselect),NhaloSeeds(:,i)
EndDo
End Program RandTest