multithreading - OpenMP block gives false results -

i appreciate point of view might did wrong using openmp. parallelized code pretty strait forward - yet single thread (i.e., call omp_set_num_threads(1)) wrong results.

i have checked intel inspector, , not have race condition, yet inspector tool indicated warning thread might approach other thread stack (i have warning in other code have, , runs openmp). not think problem.

subroutine gr(number_d, rad_d, rad_cc, spect)    use term,only: density, temperature, viscosity, water_density, &                   pressure, d_hor, d_ver, d_temp, qqq, umu   use satur,only: ff, a1, a2, aaa, bbb, sat   use delta,only: ddm, dt   use const,only: pi, g    implicit none    integer,intent(in) :: number_d   double precision,intent(in) :: rad_cc(number_d), spect(number_d)   double precision,intent(inout) :: rad_d(number_d)    double precision :: r3, dr3, c2, c0, p, q, rad_cr, sat_cr, c4, a, &                        c, d, cc, dd, cc2, dd2, rad_st, draa, dra, dm, x1                       integer ::    ddm = 0.0d0    !$omp parallel default(shared) &    !$omp private(i,r3,dr3,c2,c0,p,q,sat,sat_cr,c4,a) &   !$omp private (c,d,cc,dd,cc2,dd2,rad_st,draa,dra,dm,rad_cr,x1) &   !$omp reduction (+:ddm)   i=1,number_d      r3   = rad_cc(i)**3      dr3  = rad_d(i)**3-r3      if(dr3.lt.1.0d-100) dr3 = 1.0d-100       c2 = -dsqrt(3.0d0*bbb*r3/aaa)      c0 = -r3      p  = -0.3333333333d0*c2**2      q  = c0+0.074074074d0*c2**3      call cubic(p, q, rad_cr)      rad_cr = rad_cr - 0.3333333333d0*c2      sat_cr = dexp(aaa/rad_cr-bbb*r3/(rad_cr**3-r3))-1.0d0       dra  = dt*(sat+1.0d0-dexp(aaa/rad_drop(i)-bbb*r3/dr3))/ &            (ff*rad_d(i))       if(sat.lt.sat_cr)          if(dabs(sat).lt.1.0d-10)            p = -bbb*r3/aaa            q = -r3            call cubic(p, q, rad_st)            go 22         end if           c4  = dlog(sat+1.0d0)           = -aaa/c4         c   = (bbb-c4)*r3/c4         d   = -a*r3         p   = a*c-4.0d0*d         q   = -(a**2*d+c**2)          call cubic(p, q, x1)          cc  = dsqrt(a**2+4.d0*x1)         dd  = dsqrt(x1**2-4.d0*d)         cc2 = 0.5d0*(a-cc)         if(sat.lt.0.0d0)            dd2    = 0.5d0*(x1-dd)            rad_st = 0.5d0*(-cc2+dsqrt(cc2**2-4.0d0*dd2))         else            dd2    = 0.5d0*(x1+dd)            rad_st = 0.5d0*(-cc2-dsqrt(cc2**2-4.0d0*dd2))         end if  22       continue          draa = rad_st-rad_d(i)         if(abs(draa).lt.abs(dra))            dra = draa            dm  = 1.3333333333333333d0*pi*water_density* &                 (rad_st**3-rad_d(i)**3)         else            dm  = 4.0d0*pi*water_density*rad_d(i)**2*dra         end if         ddm = ddm+spect(i)*dm         rad_d(i) = rad_d(i) + dra      else         dm  = 4.0d0*pi*water_density*rad_d(i)**2*dra         ddm = ddm+spect(i)*dm         rad_d(i) = rad_d(i) + dra      end if    end   !$omp end parallel    return   end subroutine gr  subroutine cubic(p, q, x)    implicit none    double precision,intent(in) :: p, q   double precision,intent(out) :: x    double precision :: dis, pp, cosalfa,alfa, qq, u, v     dis = (p/3.d0)**3+(0.5d0*q)**2   if(dis.lt.0.0d0)      pp      = -p/3.0d0      cosalfa = -0.5d0*q/dsqrt(pp**3)      alfa    = dacos(cosalfa)      x       = 2.0d0*dsqrt(pp)*dcos(alfa/3.0d0)      return   else      qq  = dsqrt(dis)      u   = -0.5d0*q+qq      v   = -0.5d0*q-qq      if(u.ge.0.0d0)         u = u**0.333333333333333d0      else         u = -(-u)**0.333333333333333d0      end if      if(v.ge.0.0d0)         v = v**0.333333333333333d0      else         v = -(-v)**0.333333333333333d0      end if      x = u+v   end if  return end subroutine cubic