bug in ana_psource with openMP

[jpringle]

I am not well versed in openMP programming -- but I think I found a bug where ana_psource initializes the tracer input. I was having problems with different runs of my code producing wildly different answers from run to run under openMP. It did not do this when run in a single thread and when run under MPI. I have traced it down to the section of ana_psource.h which initializes the tracer input.

The original code (with my addition for the case JMPBUMP) went as follows

Code: Select all

!
!-----------------------------------------------------------------------
!  Set tracers point Sources and/or Sinks.
!-----------------------------------------------------------------------
!
      TRACERS : IF (ANY(LtracerSrc(:,ng))) THEN
      SOURCES(ng)%Tsrc=0.0_r8   ! initialize
!
!  Set-up tracer (tracer units) point Sources/Sinks.
!
# if defined RIVERPLUME1
        IF (DOMAIN(ng)%NorthEast_Test(tile)) THEN
          DO k=1,N(ng)
            DO is=1,Nsrc(ng)
              SOURCES(ng)%Tsrc(is,k,itemp)=T0(ng)
              SOURCES(ng)%Tsrc(is,k,isalt)=0.0_r8
            END DO
          END DO
        END IF
# elif defined JMPBUMP
        IF (DOMAIN(ng)%NorthEast_Test(tile)) THEN
          DO k=1,N(ng)
            DO is=1,Nsrc(ng)
              SOURCES(ng)%Tsrc(is,k,itemp)=T0(ng)
              SOURCES(ng)%Tsrc(is,k,isalt)=S0(ng)*0.5_r8
            END DO
          END DO
        END IF

I think the problem is with the line "SOURCES(ng)%Tsrc=0.0_r8 ! initialize". Under openMP (especially running with more tiles then threads, as is often optimal) there is no guarantee that the last thread to go through this code region is the one that initializes the sources (e.g. is the one where "(DOMAIN(ng)%NorthEast_Test(tile))" is true). Examining the tracers, sometimes they were set to zero, and sometimes they were properly set. This behavior was not deterministic.

I can get get the code to run the same way each time it is run by commenting out the line "SOURCES(ng)%Tsrc=0.0_r8 ! initialize". It serves no purpose anyway that I can tell -- you need to initialize the input tracer values to something sensible anyway, just so your model run makes sense.

Cheers, Jamie