mysterious floating point exception

[Braun__Jay]

I am working with a package provided by a university-affiliated organization.  I keep hitting an arithmetic exception on a particular line of code, and there seems to be no reason. I am wondering what actions I can take.  I am already running with -debug and no optimization. Below is the code snippet, with the "offending" line indicated in bold:

    DO k=1,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          dzp = zp(i,j,k+1)-zp(i,j,k)
          IF( runmod == 1 ) THEN
            grdscl(i,j,k) = dxdy3rd * dzp**one3rd
          ELSE IF( runmod == 2 ) THEN
            grdscl(i,j,k) = SQRT(dx*dzp)
          ELSE IF( runmod == 3 ) THEN
            grdscl(i,j,k) = SQRT(dy*dzp)
          ELSE IF( runmod == 4 ) THEN
            grdscl(i,j,k) = dzp
          END IF
        END DO
      END DO
    END DO

All variables on this line are REAL (including the 3-dimensional array).  The values are:

one3rd = 1.0/3.0 (so 0.3333333...)
dzp = -16.666626
dxdy3rd = 9.65489483

so basically the expression is taking a negative real number, raising it to a real fraction (should be fine), and the multiplying it by a positive real number.

Regarding the loop through the arrays, the indexes and the variable runmod are all correctly INTEGER.  nx is 100, ny is 400, and nz is 10.

I am running on Ubuntu 16.04.

Here is the compilation command; my only contribution was replacing -O3 wth -O0 -debug:

ifort -w -convert big_endian  -fp-model source  -I/home/braun/arps/include  -I/usr/local/hdf/include -module /home/braun/arps/modules -O0     -debug -fpe0 -check bounds -traceback -c tmix3d.f90

I'm starting to suspect a compiler option.  I tried a small test case with no options, and the result was correct.

Thanks,

Jay

[Steve_Lionel]

Um, raising a negative value to a fractional exponent results in a complex value, and hence an error for REAL values. Right?

[Braun__Jay]

Funny, I was so intent on creating a small test case, I did not consider the basic math.  My "should be fine" comment came from first getting a result on my computer's calculator and then this little program:

program j

real :: dzp, one3rd, result

one3rd = 1.0 / 3.0
dzp = -16.666626

result = dzp**one3rd
print *, result

end program j

So I compile and run it:

braun@audrey-00:~$ ifort -c j.f90
braun@audrey-00:~$ ifort -o j j.o
braun@audrey-00:~$ ./j
  -2.554363

In any case, my task is to find out why dzp is negative.  Maybe the real  "bug" is the apparently reasonable result in the small test case.

Thanks as always, Steve.

j

[Braun__Jay]

I added an additional step to the test case, multiplying that result by a real number, and got NaN.

OK, on to the next step:  where did that negative number come from?

Thanks!

j

[Braun__Jay]

Just to close this discussion -- the negative values resulted from a mismatch of parametric data.  By changing one setting that was inconsistent with the dimensions chosen for my particular data set, I eliminated the issue.  This was consistent with some documentation that I later found.

j