c-------------------------------------------------------------------- 
c                 Test OMP parallelization
c-------------------------------------------------------------------- 
      parameter (N = 1024 )
      parameter (Np = 2e7 )
      real*8 A
      COMMON /aa/ A(N+1,N+1,N+1)
      COMMON /xx/ xc(Np)
      COMMON /yy/ yc(Np)
      COMMON /zz/ zc(Np)

      real*8    sum1
      xmax = N+1.-1.e-4

      amem = (8.*N**3+3.*4.*Np)/1024.**2
      write (*,*)  ' Required Memory=',amem,'Mb'
 		
C$OMP PARALLEL DO DEFAULT(SHARED) 
C$OMP+PRIVATE (i,j,k)
      Do k=1,N+1
      Do j=1,N+1
         Do i=1,N+1
            A(i,j,k) = 0.
         EndDo 
      EndDo 
      EndDo 
      Nseed = 121091
      Do ip=1,Np
         xc(ip) =1.+RANDd(Nseed)*n
         xc(ip) =MIN(xc(ip),xmax)
         If(xc(ip).ge.N+1)write(*,*) ' error:',xc(ip),xx,Nseed
      EndDo       
      Do ip=1,Np
         yc(ip) = 1.+RANDd(Nseed)*N
         yc(ip) =MIN(yc(ip),xmax)
      EndDo       
      Do ip=1,Np
         zc(ip) = 1.+RANDd(Nseed)*N
         zc(ip) =MIN(zc(ip),xmax)
      EndDo

      Do iloop =1,1000
C$OMP PARALLEL DO DEFAULT(SHARED) 
C$OMP+PRIVATE (ip)
      Do ip=1,Np
         CALL AddDens(xc(ip),xc(ip),xc(ip))
      EndDo 
      EndDo

      sum1 =0.
C$OMP PARALLEL DO DEFAULT(SHARED) 
C$OMP+PRIVATE (i,j,k)
C$OMP+REDUCTION(+:sum1)
      Do k=1,N+1
      Do j=1,N+1
         Do i=1,N+1
            sum1 =sum1 + A(i,j,k)**2
         EndDo 
      EndDo 
      EndDo 

      write(*,*) ' Np=', Np,' N=',N,' Sum=',sum1

      stop
      end
!--------------------------------------------------------------
      SUBROUTINE AddDens(x,y,z)
      parameter (N = 1024 )
      real*8 A
      COMMON /aa/ A(N+1,N+1,N+1)

      i = int(x) 
      j=  int(y)
      k= int(z)
      dx = x-i
      dy = y-j
      dz = z-k
      gx = 1.-dx
      gy = 1.-dy
      gz = 1.-dz
C$OMP ATOMIC
      A(i,j,k) = A(i,j,k) +gx*gy*gz
C$OMP ATOMIC
      A(i+1,j,k) = A(i+1,j,k) +dx*gy*gz
C$OMP ATOMIC
      A(i,j+1,k) = A(i,j+1,k) +gx*dy*gz
C$OMP ATOMIC
      A(i+1,j+1,k) = A(i+1,j+1,k) +dx*dy*gz
C$OMP ATOMIC
      A(i,j,k+1) = A(i,j,k+1) +gx*gy*dz
C$OMP ATOMIC
      A(i+1,j,k+1) = A(i+1,j,k+1) +dx*gy*dz
C$OMP ATOMIC
      A(i,j+1,k+1) = A(i,j+1,k+1) +gx*dy*dz
C$OMP ATOMIC
      A(i+1,j+1,k+1) = A(i+1,j+1,k+1) +dx*dy*dz

      return
      end

C------------------------------------------------
C				                                       random number generator
      FUNCTION RANDd(M)
C------------------------------------------------
      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