Is There a Better Double-Precision Assignment in Fortran 90?

In Fortran 90 (using gfortran on Mac OS X) if I assign a value to a double-precision variable without explicitly tacking on a kind, the precision doesn't "take." What I mean is, if I run the following program:

program sample_dp

implicit none

integer, parameter :: sp = kind(1.0)
integer, parameter :: dp = kind(1.0d0)

real(sp) :: a = 0.
real(dp) :: b = 0., c = 0., d = 0.0_dp, e = 0_dp

! assign values
a = 0.12345678901234567890
b = 0.12345678901234567890
c = DBLE(0.12345678901234567890)
d = 0.12345678901234567890_dp

write(*,101) a, b, c, d
101 format(1x, 'Single precision: ',  T27, F17.15, / &
           1x, 'Double precisison: ', T27, F17.15, / &
           1x, 'Double precision (DBLE): ', T27, F17.15, / &
           1x, 'Double precision (_dp): ',  T27, F17.15)

end program

I get the result:

Single precision:        0.123456791043282
Double precision:        0.123456791043282
Double precision (DBLE): 0.123456791043282
Double precision (_dp):  0.123456789012346

The single precision result starts rounding off at the 8th decimal as expected, but only the double precision variable I assigned explicitly with _dp keeps all 16 digits of precision. This seems odd, as I would expect (I'm relatively new to Fortran) that a double precision variable would automatically be double-precision. Is there a better way to assign double precision variables, or do I have to explicitly type them as above?

Solution 1:

A real which isn't marked as double precision will be assumed to be single precision. Just because sometime later you assign it to a double precision variable, or convert it to double precision, that doesn't mean that the value will 'magically' be double precision. It doesn't look ahead to see how the value will be used.

Solution 2:

There are several questions linking here so it is good to state some details more explicitly with examples, especially for beginners.

As stated by MRAB in his correct answer, an expression is always evaluated without any context, so

 0.12345678901234567890

is a default (single) precision floating literal, no matter where does it appear. The same holds to floating point numbers in the exponential form

 0.12345678901234567890E0

it is also a default precision number.

If one want to use a double precision constant, one can use D instead of E in the above form. Even if such a double precision constant is assigned to a default precision variable, it is first treated as a double precision number and then it is converted to default precision.

The way you are using in your question (employing the kind notation and several kind constants) is more general and more modern, but the principle is the same.

 0.12345678901234567890_sp

is a number of kind sp and

 0.12345678901234567890_dp

is a number of kind dp and it does not matter where do they appear.

As your example shows, it is not only about assignment. In the line

 c = DBLE(0.12345678901234567890)

first the number 0.12345678901234567890 is default precision. Then it is converted to double precision by DBLE, but that is done after some of the digits are already lost. Then this new double precision number is assigned to c.