SIGFPE error with gfortran 4.8.5 handling

Question

I am using a computational fluid dynamics software that is compiled with gfortran version 4.8.5 on Ubuntu 16.04 LTS. The software can be compiled with either single precision or double precision and the -O3 optimization option. As I do not have the necessary computational resources to run the CFD software on double precision I am compiling it with single precision and the following options

ffpe-trap=invalid,zero,overflow

I am getting a SIGFPE error on a line of code that contains the asin function-

 INTEGER, PARAMETER :: sp = SELECTED_REAL_KIND( 6, 37) !< single precision
 INTEGER, PARAMETER :: wp = sp
 REAL(KIND=wp) zsm(:,:)

ela(i,j) = ASIN(zsm(ip,jp))

In other words the inverse sin function and this code is part of a doubly nested FOR loop with jp and ip as the indices. Currently the software staff is unable to help me for various other reasons and so I am trying to debug this on my own. The SIGFPE error is only being observed in the single precision compilation not double precision compilation.

I have inserted the following print statements in my code prior to the line of code that is failing i.e. the asin function call. Would this help me with unraveling the problem that I am facing ? This piece of code is executed for every time step and it is occurring after a series of time steps. Alternatively what other steps can I do to help me fix this problem ? Would adding "precision" to the compiler flag help ?

  if (zsm(ip,jp) >= 1.0 .or. zsm(ip,jp) <= -1.0) then
       print *,zsm(ip,jp),ip,jp
  end if

EDIT I took a look at this answer Unexpected behavior of asin in R and I am wondering whether I could do something similar in fortran i.e. by using the max function. If it goes below -1 or greater than 1 then round it off in the proper manner. How can I do it with gfortran using the max function ?

On my desktop the following program executes with no problems(i.e. it has the ability to handle signed zeros properly) and so I am guessing the SIGFPE error occurs with either the argument greater than 1 or less than -1.

 program testa

 real a,x

 x = -0.0000
 a = asin(x)
 print *,a
 end program testa

Answer 1

We have min and max functions in Fortran, so I think we can use the same method as in the linked page, i.e., asin( max(-1.0,min(1.0,x) ). I have tried the following test with gfortran-4.8 & 7.1:

program main
    implicit none
    integer, parameter :: sp = selected_real_kind( 6, 37 )
    integer, parameter :: wp = sp
    ! integer, parameter :: wp = kind( 0.0 )
    ! integer, parameter :: wp = kind( 0.0d0 )
    real(wp) :: x, a

    print *, "Input x"
    read(*,*) x

    print *, "x =", x
    print *, "equal to 1 ? :", x == 1.0_wp
    print *, asin( x )
    print *, asin( max( -1.0_wp, min( 1.0_wp, x ) ) )
end

which gives with wp = sp (or wp = kind(0.0) on my computer)

$ ./a.out
 Input x
1.00000001
 x =   1.00000000    
 equal to 1 ? : T
   1.57079625             (<- 1.5707964 for gfortran-4.8)
   1.57079625    

$ ./a.out
 Input x
1.0000001
 x =   1.00000012    
 equal to 1 ? : F
              NaN
   1.57079625 

and with wp = kind(0.0d0)

$ ./a.out
 Input x
1.0000000000000001
 x =   1.0000000000000000     
 equal to 1 ? : T
   1.5707963267948966
   1.5707963267948966     

$ ./a.out
 Input x
1.000000000000001     
 x =   1.0000000000000011     
 equal to 1 ? : F
                       NaN
   1.5707963267948966 

---

If it is necessary to modify a lot of asin(x) and the program relies on a C or Fortran preprocessor, it may be convenient to define some macro like

#define clamp(x) max(-1.0_wp,min(1.0_wp,x))

and use it as asin( clamp(x) ). If we want to remove such a modification, we can simply change the definition of clamp() as #define clamp(x) (x). Another approach may be to define some asin2(x) function that limits x to [-1,1] and replace the built-in asin by asin2 (either as a macro or a Fortran function).