I have a parallelised MPI code which solves non-linear PDEs (basically LES simulations), written in Fortran 90. In a particular section of the code I have to multiply double precision variables with an other double precision variable.

u(iwrap_p(:),j,k) = alpha(:) * u(iwrap_p(:),j,k) + beta(:) * u(iwrap_p(:),j,k)

where alpha and beta are double precision variables whose values vary between 0. and 1.0, u is the velocity,and the sum of alpha and beta at any position is 1.0. The output of this operation should be equal to u(iwrap_p(:), j, k). But, performing this operation i.e. multiplication of the velocity with alpha and beta introduces round-off errors which keep on growing with time and finally the solution becomes non-physical.

Alpha and beta are calculated by the following,

beta = 0.5_rprec * ( 1.0_rprec - cos (pi * real (i - istart, rprec) & / real(iend - istart, rprec)) )

alpha = 1._rprec - beta

where 'rprec' is the user defined double precision datatype.

I am using the following compiler arguments

FFLAGS = -O3 -fp-model source

I have tried compiling with and without optimisations, but nothing seem to work.

The error disappears when alpha and beta are declared as single precision variables. What might be the cause of this error?