OpenMP inconsistently segfaulting

I parallelized my code with openMP, and now have this bug in my code that’s really odd. If it’s the first time I run on a computer, it segfaults, but if I run the program again it runs fine. I’ve never seen anything like this and was curious if anyone has seen this before or has good resources or hints on how to fix this. It’s a very simple loop that it segfaults on, and I parallelized every loop with default none. The variable declaration followed by the loop is below. I had previously posted a stripped down version of the loop for readability, but it was requested that I fill it in for the sake of more clarity :

 real(dp)    :: q, oneOverq, L, D, Pressure, u, v, En, rho, oneOverrho
 real(dp)    :: normal(Mesh%numDims), wL(Mesh%numFields), cd_dql(Mesh%numFields, Mesh%numEdges)
 real(dp)    :: norm_dx1(Mesh%numDims, Mesh%numDims), norm_dx2(Mesh%numDims, Mesh%numDims)
 real(dp)    :: L_dx(Mesh%numDims, Mesh%numNodes), D_dx(Mesh%numDims, Mesh%numNodes)
 real(dp)    :: p_dql(Mesh%numFields), L_dql(Mesh%numFields, Mesh%numEdges), Res(Mesh%numFields, Mesh%numTri)
 INTEGER(i4) :: i, j, k, m, nodecount, cellNum, edgeNum, Node1, Node2, node, cell, cells(3)


   !$OMP PARALLEL DO Reduction(+:L, D, cL_dql, L_dx, cD_dql, D_Dx) default(none) shared(Mesh, W, grad, gamma, dRdX, Jac_2nd) &
   !$OMP private(cellnum, edgenum, node1, node2, rho, oneoverrho, u, v, En, pressure, normal, norm_dx1, norm_dx2, p_dql, wL, cell, node, cells)
   !DIR$ IVDEP:LOOP
   do cellnum = 1, Mesh%numTri
       edgeNum = Mesh%airfoilEdges(i)
       Node1   = Mesh%edgelist(edgenum,1)
       Node2   = Mesh%edgelist(edgenum,2)
       cellNum = Mesh%edgeList(edgenum,3)
       call cell_stencil(Mesh, cellnum, cells)

       wL = W(:, cellnum)

       call linear_limited_reconstruction(cellnum, edgenum, Mesh, grad(:,:,cellNum), wL)


       rho = wL(1)
       oneOverrho = 1.0_dp/rho
       u   = wL(2)*oneOverrho
       v   = wL(3)*oneOverrho
       En  = wL(4)*oneOverrho

       call get_pressure(rho, En, u, v, pressure)

       p_dql(1) = (gamma - 1._dp)*half*( u**2 + v**2 )
       p_dql(2) = (gamma - 1._dp)*-u
       p_dql(3) = (gamma - 1._dp)*-v
       p_dql(4) = gamma-1._dp


       call norm_dx(node1, node2, Mesh, normal, norm_dx1, norm_dx2)

       L     = L + Pressure*normal(2)
       D     = D + Pressure*normal(1)


       do j = 1, 4
           if (j == 1) then 
               cell = cellnum
           else
               cell = cells(j-1)
           end if
           if (cell /= 0) then
               cL_dql(:,cell) = cL_dql(:,cell) + matmul(P_dql(:),Jac_2nd%wL_rc_dq(:,:,j,edgenum))*normal(2)
               cD_dql(:,cell) = cD_dql(:,cell) + matmul(P_dql(:),Jac_2nd%wL_rc_dq(:,:,j,edgenum))*normal(1)
           end if 
       end do

       L_dx(:, node1)   = L_dx(:, node1) + Pressure*norm_dx1(2, :)
       L_dx(:, node2)   = L_dx(:, node2) + Pressure*norm_dx2(2, :)
       D_dx(:, node1)   = D_dx(:, node1) + Pressure*norm_dx1(1, :)
       D_dx(:, node2)   = D_dx(:, node2) + Pressure*norm_dx2(1, :)

       do j = 1, 6
           node = Mesh%NodeStencil(j, cellnum)
           if (node /= 0) then
               do k = 1, Mesh%numFields 
                   L_dx(:,node) = L_dx(:,node) + P_dql(k)*dRdX%wL_rc_dx(k,:,j,edgenum)*normal(2)              
                   D_dx(:,node) = D_dx(:,node) + P_dql(k)*dRdX%wL_rc_dx(k,:,j,edgenum)*normal(1)
               end do
           end if
       end do  

   end do
   !$OMP END PARALLEL DO

I can put in the full code up in the code block, but for the sake of readability I thought it was best not to. If anyone has any advice I’d greatly appreciate it.

Thanks!