openmp并行Fortran后每次结果不一致

[? Vita]

进行并行化之后每次结果不一致，对赋值代码进行了atomic保护也不好使。有没有大佬帮忙看看怎么回事。
代码很长。也不知道有没有大佬有耐心看完....

  implicit none
  integer n,m,node,node1,node2,ni,mj,p,q,k
  integer i,j

  real(4) rcm(3),rcn(3),dist
  real r1f(3),r2f(3),r3f(3),nSm(3),rjf(3,SNum_samp),r1s(3),r2s(3),r3s(3),ris(3,SNum_samp)
  real Am,An,lj,li,rmj(3),rni(3),rtmpv(3)
  real Jrm(3,maxRWG),Irn(3,maxRWG),Jlengthm(maxRWG),Ilengthn(maxRWG)

  complex coef1,coef2
  complex Amn_E,Zji,I1ss,I2ss(3),I3ss(3),I4ss
  complex ctmp,ctmp1
    
  real,external:: re_Norm
  integer*4 time1,time2,time
  
  call SYSTEM_CLOCK(time1)
  write(*,*) "Start to compute the SS part of the impedance matrix, on the way..."
  
  allocate(SparseValue(Num_side,Num_side))
  SparseValue=0
      
  coef1=cmplx(0.0,omega*MU0/(16.0*PI))
  coef2=cmplx(0.0,1.0/(4.0*PI*omega*EPS0))
 
!$OMP parallel do default(none) shared(Num_Triangle,triangle_centroid,triangle_area,point_cor,triangle_point,triangle_rwgs,side_index,side_length,SparseValue,threshold_R,coef1,coef2) private(n,m,node,node1,node2,ni,mj,i,j,rcm,rcn,dist,r1f,r2f,r3f,rjf,r1s,r2s,r3s,ris,Am,An,lj,li,rmj,rni,Jrm,Irn,Jlengthm,Ilengthn,Amn_E,ctmp,I1ss,I2ss,I3ss,I4ss)
  do m=1,Num_Triangle    
    rcm=triangle_centroid(:,m)
	   Am=triangle_area(m)
 
    r1f=point_cor(:,triangle_point(1,m))
    r2f=point_cor(:,triangle_point(2,m))
    r3f=point_cor(:,triangle_point(3,m))

	   do k=1,SNum_samp
      rjf(:,k)=epsilS(k)*r1f+etaS(k)*r2f+xiS(k)*r3f
    enddo
 
    do mj=1,triangle_rwgs(m)%N_rwg
      j=triangle_rwgs(m)%rwg(mj)
        
      if(side_index(3,j)==m) then
        !flagm=1.0
        Jlengthm(mj)=side_length(j)
      else  !if(side_index(4,j)==m) then
        !flagm=-1.0
        Jlengthm(mj)=-side_length(j)
      endif
        
   	  node1=side_index(1,j); node2=side_index(2,j)
      do p=1,3
        node=triangle_point(p,m)
        if((node/=node1).and.(node/=node2)) exit
      enddo         
      Jrm(:,mj)=point_cor(:,node)   	          
    enddo 
 
    do n=1,Num_Triangle
      rcn=triangle_centroid(:,n)
      An=triangle_area(n)

      r1s=point_cor(:,triangle_point(1,n))     !source triangle coordinates
      r2s=point_cor(:,triangle_point(2,n))
      r3s=point_cor(:,triangle_point(3,n))

      do k=1,SNum_samp
	       ris(:,k)=epsilS(k)*r1s+etaS(k)*r2s+xiS(k)*r3s
      enddo
                    
      dist=re_Norm(rcm-rcn)          
      
      call Comp_eltsI(I1ss,I2ss,I3ss,I4ss,dist>=threshold_R,An,rjf,ris,r1s,r2s,r3s,SNum_samp,SNum_samp,weightS,weightS)
       
      do ni=1,triangle_rwgs(n)%N_rwg
        i=triangle_rwgs(n)%rwg(ni)
		      if(side_index(3,i)==n) then
	         !flagn=1.0
          Ilengthn(ni)=side_length(i)
        else  !if(side_index(4,i)==n) then
	         !flagn=-1.0
          Ilengthn(ni)=-side_length(i)
	       endif
              
	       node1=side_index(1,i); node2=side_index(2,i)
	       do q=1,3
          node=triangle_point(q,n)
          if((node/=node1).and.(node/=node2)) exit
	       enddo
        Irn(:,ni)=point_cor(:,node)
      enddo  
      
      do mj=1,triangle_rwgs(m)%N_rwg
        j=triangle_rwgs(m)%rwg(mj)
        rmj=Jrm(:,mj)       
        lj=Jlengthm(mj) 
 
        ctmp=coef1*(I4ss-dot_product(rmj,I3ss))-coef2*I1ss
         
        do ni=1,triangle_rwgs(n)%N_rwg
          i=triangle_rwgs(n)%rwg(ni)
                    
          rni=Irn(:,ni)               
          li=Ilengthn(ni)
                    
          Amn_E=lj*li*(coef1*(dot_product(rmj,rni)*I1ss-dot_product(rni,I2ss))+ctmp)

!$OMP atomic
          SparseValue(i,j)=SparseValue(i,j)+Amn_E     
         
	       enddo
      enddo
    enddo

	 enddo
!$OMP end parallel do

[境�]

请问解决了吗，我也是这样的问题