I don't have enough rep to comment so I'll "answer" regarding Rick Thompson's suggestion for testing infinity.
if (A-1 .eq. A)
This will also be true if A is a very large floating point number, and 1
is below the precision of A.
A simple test:
subroutine test_inf_1(A)
real, intent(in) :: A
print*, "Test (A-1 == A)"
if (A-1 .eq. A) then
print*, " INFINITY!!!"
else
print*, " NOT infinite"
endif
end subroutine
subroutine test_inf_2(A)
real, intent(in) :: A
print*, "Test (A > HUGE(A))"
if (A > HUGE(A)) then
print*, " INFINITY!!!"
else
print*, " NOT infinite"
endif
end subroutine
program test
real :: A,B
A=10
print*, "A = ",A
call test_inf_1(A)
call test_inf_2(A)
print*, ""
A=1e20
print*, "A = ",A
call test_inf_1(A)
call test_inf_2(A)
print*, ""
B=0.0 ! B is necessary to trick gfortran into compiling this
A=1/B
print*, "A = ",A
call test_inf_1(A)
call test_inf_2(A)
print*, ""
end program test
outputs:
A = 10.0000000
Test (A-1 == A)
NOT infinite
Test (A > HUGE(A))
NOT infinite
A = 1.00000002E+20
Test (A-1 == A)
INFINITY!!!
Test (A > HUGE(A))
NOT infinite
A = Infinity
Test (A-1 == A)
INFINITY!!!
Test (A > HUGE(A))
INFINITY!!!