Make constants explicitly double precision, otherwise they are unintentionally treated as single precision by Fortran 90

[n8xm]

When making an assignment statement such as the below, one might think that pi is treated as a double precision float:

real(kind=dp),public,parameter :: pi=3.141592654

However, Fortran 90, unlike Fortran 77, will treat pi as a single precision float, rather than a double. This is a known (and rather counterintuitive) behavior of Fortran 90. Below is a minimum working example of this behavior in Fortran:

program flt_vs_double                                                                                                                                                        
        use, intrinsic :: iso_fortran_env, only: sp=>real32, dp=>real64                                                                                    
        implicit none                                                                                                                                                  
                                                                                                                                                                                      
        real(dp) :: foo                                                                                                                                                              
        real(dp) :: bar                                                                                                                                                                                  
        foo = 3.141592654_dp                                                                                                                                                                             
        bar = 3.141592654                                                                                                                                   
                                                                                                                                                          
        print *, 'foo^3 is:', foo*foo*foo                                                                                                                 
        print *, 'bar^3 is:', bar*bar*bar                                                                                                                 
                                                                                                                                                          
end program flt_vs_double    

The output of the above program when I run it (compiling with gfortran, as part of GCC 11.1.0) is:

 foo^3 is:   31.006276692445560     
 bar^3 is:   31.006279268784656

Note the difference in results.

Now, consider the following C program:

#include <stdio.h>   
                                                                       
void main(){         
  double foo = 3.141592654;
  float bar = 3.141592654;
  printf(" foo^3 is:   %.15f\n", foo*foo*foo);
  printf(" bar^3 is:   %.15f\n", bar*bar*bar);
} 

The output of the above is:

 foo^3 is:   31.006276692445560
 bar^3 is:   31.006278991699219

Notice how foo^3 matches the Fortran results, and bar^3 matches the Fortran results up to 6 decimal places. This makes sense given single precision (a rule of thumb for single precision is 7 decimal places when using scientific notation). In other words, these results demonstrate that Fortran treats bar as a single precision float, because we did not append the _dp to the end of it.

See also the following sources online:

• https://stackoverflow.com/a/6146102
• https://www.fortran90.org/src/best-practices.html#floating-point-numbers
• https://www.fortran90.org/src/gotchas.html#floating-point-numbers-gotcha
• https://stackoverflow.com/a/45156962

In this pull request, I have appended _dp to constants used by this code, in places where it appears that the intention was to use double precision. Of course, this subtly changes the numerical results. I have tested with test1.in and it does appear to work.

Pinging contributors: @garrelt @abken601