c-------------------------------------------------------------------- 
c                 Test OMP parallelization
c-------------------------------------------------------------------- 
      parameter (N = 1200 )
      COMMON /aa/ A(N,N,N)
      COMMON /bb/ B(N,N,N)
      real*8    sum,sum1

      write (*,*)  ' Required Memory=',8.*float(N)**3/1024.**2,'Mb'

C$OMP PARALLEL DO DEFAULT(SHARED) 
C$OMP+PRIVATE (i,j,k)
      Do j=1,N
      Do i=1,N
         Do k=1,N
            A(i,k,j)=0
         EndDo 
      EndDo 
      EndDo 

C$OMP PARALLEL DO DEFAULT(SHARED) 
C$OMP+PRIVATE (i,j,k)
      Do j=1,N
      Do i=1,N
         Do k=1,N
            B(i,k,j)=0
         EndDo 
      EndDo 
      EndDo 

      sum = 0.
C$OMP PARALLEL DO DEFAULT(SHARED) 
C$OMP+PRIVATE (i,j,k)
C$OMP+REDUCTION(+:sum)
      Do j=1,N
      Do i=1,N
         Do k=1,N
            A(i,k,j) = ((mod(i,5)+mod(j,5)+mod(k,5))/100.)
            sum =sum + A(i,k,j)**2
         EndDo 
      EndDo 
      EndDo 
      write(*,*) ' Sum of elements of Matrix a = ',sum
      m=3
      CALL SMOOTH(m)
      sum1 = 0.
C$OMP PARALLEL DO DEFAULT(SHARED) 
C$OMP+PRIVATE (i,j,k)
C$OMP+REDUCTION(+:sum1)
      Do j=1,N
      Do i=1,N
         Do k=1,N
            sum1 =sum1 + A(i,k,j)**2
         EndDo 
      EndDo 
      EndDo 
      write (*,*)  ' smoothing=',m
      write (*,*)  ' average=',sqrt(sum/N**2)
      write (*,*)  ' average=',sqrt(sum1/N**2)
           stop
           end

c-------------------------------------------------------------------- 
c
      SUBROUTINE SMOOTH(m)
      parameter (N = 1200 )
      COMMON /aa/ A(N,N,N)
      COMMON /bb/ B(N,N,N)
      real*8    sum,norm

      norm =(2.*m+1)**3
      write (*,*)  ' norm=',norm,m
C$OMP PARALLEL DO DEFAULT(SHARED) 
C$OMP+PRIVATE (i,j,k,sum,i1,j1,k1,ii,jj,kk)
      Do k=1,N,2
         Do j=1,N,2
            Do i=1,N,2
               sum =0.
               Do k1=-m,m
               Do j1=-m,m
               Do i1=-m,m
                 ii = min(max(i+i1,1),N)
                 jj = min(max(j+j1,1),N)
                 kk = min(max(k+k1,1),N)
                 sum = sum+A(ii,jj,kk)
               EndDo 
               EndDo 
               EndDo 
               B(i,j,k) =sum/norm
            EndDo 
         EndDo 
      EndDo 
      write(*,*) ' end first loop in smooth'
C$OMP PARALLEL DO DEFAULT(SHARED) 
C$OMP+PRIVATE (i,j,k)
      Do j=1,N
      Do i=1,N
         Do k=1,N
            A(i,j,k)=B(i,j,k)
         EndDo 
      EndDo 
      EndDo 
      write(*,*) ' done smooth'
      Return
      end