Q: OpenMP loop with if-statement only causing NaN I found some NaN's in a data array in my Fortran code, and may have isolated the issue to an OpenMP do loop. When this loop runs, NaNs result in sp%ptl(i)%ph(6):
!$OMP PARALLEL DO &
!$OMP PRIVATE( ITH, I )
do ith=1, sml_nthreads
  do i=i_beg(ith), i_end(ith)
    if(sp%ptl(i)%ph(3) >= 2pi .or. sp%ptl(i)%ph(3)< 0D0 ) then
      sp%ptl(i)%ph(3) = modulo(sp%ptl(i)%ph(3),2pi)
    endif
  enddo
enddo

But if I run the same, but add a line adding 0D0 to a dummy variable, the NaN's go away in sp%ptl(i)%ph(6):
!$OMP PARALLEL DO &
!$OMP PRIVATE( ITH, I )
do ith=1, sml_nthreads
  do i=i_beg(ith), i_end(ith)
    if(sp%ptl(i)%ph(3) >= 2pi .or. sp%ptl(i)%ph(3)< 0D0 ) then
      sp%ptl(i)%ph(3) = modulo(sp%ptl(i)%ph(3),2pi)
    endif

    tmp = tmp + 0D0
  enddo
enddo

Of course, there is more to the real code and this is not a minimal working example. My question though is why does adding any line to the do loop cause that the sp%ptl(i)%ph(6) never get the NaN's? Is having only an if-statement in the OpenMP inner do loop a bad idea? It's a confusing issue to me right now how this is working.

UPDATE Here is minimal example, not quite working as the larger code yet, as it doesn't have NaN's but rather large numbers in random points of the array ptl, but at the very least shows the basic workflow. I compile the same as with the larger codebase (Intel compiler, 18.0.1.163), then run this with srun -n1 -c24 . 
I did some further tests with the larger codebase, and found that recompiling the subroutines represented here by "mymod" without optimization (i.e. -O0 -g -C) makes the NaN's go away. 
UPDATE 2 Never mind about the large numbers, I had simply forgotten the initialization of ptl(i)%ph (now added at beginning of push), now that its added I never get NaN or  large numbers in this minimal example (still there in the larger code with optimization on). 
module mymod
    integer, parameter :: ptl_nphase=8
    integer, parameter :: num=1000
    integer, parameter :: sml_nthreads=24
    type ptl_type
        real(8) :: ph(ptl_nphase)
    end type ptl_type
contains
    logical function is_nan(a)
        implicit none
        real (8) :: a

        is_nan = .not. ( a > 1D0 .or. a < 2D0 )

    end function is_nan

    subroutine split_indices(total,num_pieces,ibeg,iend)
        implicit none

        integer :: total
        integer :: num_pieces
        integer :: ibeg(num_pieces), iend(num_pieces)
        integer :: itmp1, itmp2, ioffset, i

        if (num_pieces > 0) then
            itmp1 = total/num_pieces
            itmp2 = mod(total,num_pieces)
            ioffset = 0
            do i=1,itmp2
            ibeg(i) = ioffset + 1
            iend(i) = ioffset + (itmp1+1)
            ioffset = iend(i)
            enddo
            do i=itmp2+1,num_pieces
            ibeg(i) = ioffset + 1
            if (ibeg(i) > total) then
                iend(i) = ibeg(i) - 1
            else
                iend(i) = ioffset + itmp1
                ioffset = iend(i)
            endif
            enddo
        endif
    end subroutine split_indices


    subroutine calc_source(ptl,icycle)
        implicit none
        type(ptl_type) :: ptl(num)
        integer :: ith, i, i_beg(sml_nthreads), i_end(sml_nthreads)
        integer :: icycle

        call split_indices(num, sml_nthreads, i_beg, i_end)
        !$OMP PARALLEL DO &
        !$OMP PRIVATE( ITH, I )
        do ith=1, sml_nthreads
        do i=i_beg(ith), i_end(ith)
        ptl(i)%ph(6) = ptl(i)%ph(6) + 1D0
        enddo
        enddo

        if (icycle==1) then
            !$OMP PARALLEL DO &
            !$OMP PRIVATE( ITH, I )
            do ith=1, sml_nthreads
                do i=i_beg(ith), i_end(ith)
                    ptl(i)%ph(7) = ptl(i)%ph(6)
                enddo
            enddo
        endif

    end subroutine calc_source


    subroutine push1(ptl)
        implicit none
        type(ptl_type) :: ptl(num)
        integer :: ith, i, i_beg(sml_nthreads), i_end(sml_nthreads)
        real(8) :: arr1(5)

        call split_indices(num, sml_nthreads, i_beg, i_end)
        !$OMP PARALLEL DO &
        !$OMP PRIVATE( ITH, I )
        do ith=1, sml_nthreads
            do i=i_beg(ith), i_end(ith)
                call random_number(arr1)
                ptl(i)%ph(1:5) = ptl(i)%ph(1:5) + arr1
            enddo
        enddo

    end subroutine push1


    subroutine push(ptl)
        implicit none
        type(ptl_type) :: ptl(num)
        integer :: icycle
        integer :: ith, i, i_beg(sml_nthreads), i_end(sml_nthreads)

        call split_indices(num, sml_nthreads, i_beg, i_end)
        !$OMP PARALLEL DO &
        !$OMP PRIVATE( ITH, I )
        do ith=1, sml_nthreads
            do i=i_beg(ith), i_end(ith)
               ptl(i)%ph(:) = 0D0
            enddo
        enddo

        do icycle=1,100
            call calc_source(ptl,icycle)

            call push1(ptl)

            call split_indices(num, sml_nthreads, i_beg, i_end)
            !$OMP PARALLEL DO &
            !$OMP PRIVATE( ITH, I )
            do ith=1, sml_nthreads
                do i=i_beg(ith), i_end(ith)
                    ptl(i)%ph(3) = modulo(ptl(i)%ph(3),6.28)
                enddo
            enddo

        enddo

    end subroutine push

end module mymod


program main
    use mymod
    implicit none
    type(ptl_type) :: ptl(num)
    integer :: ith, i, i_beg(sml_nthreads), i_end(sml_nthreads)

    call push(ptl)

    !check for nan
    call split_indices(num, sml_nthreads, i_beg, i_end)
    !$OMP PARALLEL DO &
    !$OMP PRIVATE( ITH, I )
    do ith=1, sml_nthreads
        do i=i_beg(ith), i_end(ith)
!            if (is_nan(ptl(i)%ph(6))) then
                !print *,'is_nan',i
                print *,ptl(i)%ph(6)
!            endif
        enddo
    enddo

end program main


A: I found a solution, but still without understanding why, and still without being able to generate a minimal, working example (my simple example never replicated the issue in the larger code). This has to do with the optimization, and how loops are inlined and vectorized, but when looking through the optimization reports, there wasn't a clear difference when compiling with the tmp = tmp+0D0 line and not. 
So I decided to start removing code. If I removed the call to calc_source, the NaN's went away. When I put it back, and just used a blank calc_source subroutine, NaNs came back. I moved calc_source to be in the same file as push, and the NaN's went away.
