PROGRAM REPRODUCER

    ! Build with: ifort  -O2 -fp-model fast=1 -xAVX -g -qopenmp t.f90

    USE OMP_LIB
    IMPLICIT NONE

    INTEGER, PARAMETER        :: MHNOD = 100
    REAL   , DIMENSION(MHNOD) :: A, SS, ALFA, ROF, ROG, VV0

    INTEGER                   :: K
    REAL, DIMENSION(MHNOD)    :: CAUX, CAUXT, DAUX, DAUXT
    REAL                      :: TERM1, TERM2

    ! Initialize
    DO K=1,MHNOD
       TERM1   = 0.001 * K 
       ALFA(K) = 0.5 + TERM1
       A(K)    = 0.2
       SS(K)   = 1.2 + TERM1
       ROF(K)  = 0.7 - TERM1
       ROG(K)  = ROF(K) * 0.001
       VV0(K)  = 1.0
    ENDDO

    ! Run 8 threads
    call OMP_SET_NUM_THREADS(8)

!$omp parallel do private(K,TERM1,TERM2)
    DO K=2,78
       TERM1=SS(K)*ALFA(K)/(1.0+ALFA(K)*(SS(K)-1.0))
       TERM2=ROF(K)-ROG(K)
       CAUXT(K)=(ROF(K)-TERM1*TERM2)/A(K)
       DAUXT(K)=VV0(K)*ALFA(K)*TERM2*(1.0-TERM1)
    ENDDO

    ! Run 1 thread
    call OMP_SET_NUM_THREADS(1)

!$omp parallel do private(K,TERM1,TERM2)
    DO K=2,78
       TERM1=SS(K)*ALFA(K)/(1.0+ALFA(K)*(SS(K)-1.0))
       TERM2=ROF(K)-ROG(K)
       CAUX(K)=(ROF(K)-TERM1*TERM2)/A(K)
       DAUX(K)=VV0(K)*ALFA(K)*TERM2*(1.0-TERM1)
    ENDDO

    ! Check results
    DO K=2,78
       print '(a,i6,2f14.8)','CAUX: ', K, CAUX(K), CAUXT(K)
       IF (CAUX(K) /= CAUXT(K)) THEN
          print *,'CAUX DIFF: ', K
       ENDIF
    ENDDO

END PROGRAM
