PRACE NPT

      Module ran_module

! ----------------------------------------------------------------------

      Use                    numerics

      Integer, Parameter  :: m1 = 259200, a1 = 7141, c1 = 54773,

     &                       m2 = 134456, a2 = 8121, c2 = 28411

      Real(l_), Parameter :: one = 1.0_l_, rm1 = one/m1,

     &                       rm2 = one/m2

      Integer             :: x1, x2

! ----------------------------------------------------------------------

      Contains

      Function dran0()                               Result( ran )

! -----------------------------------------------------------------

      Implicit   None

! -----------------------------------------------------------------

! --- dran0 returns a uniform deviate in [0,1).

!

! --- The algorithm is loosely based on an algorithm from

!     Press & Teukolsky et.al. and based on the linear congruential

!     method with choices for M, A, and C that are given by

!     D. Knuth in "Semi-numerical algorithms".

! --- Input parameters:

!     Integer  - a1, c1, m1, a2, c2, m2. The parameters of the two

!                linear congruent relations used. They are passed

!                via module 'ran_module'.

!     Integer  - x1, x2. Seeds for the two linear congruences.

!                Passed via module 'ran_module'.

!

! --- Output-parameters:

!     Real(l_) - ran.  Uniform deviate in [0,1)

! ------------------------------------------------------------------

!

      Real(l_)            :: ran

! ------------------------------------------------------------------

      x1  = Mod( a1*x1 + c1, m1 )

      x2  = Mod( a2*x2 + c2, m2 )

      ran = ( Real( x1, l_ ) + Real( x2, l_ )*rm2 )*rm1

! -----------------------------------------------------------------

      End Function dran0

      Subroutine ranfil( a, n )

! ----------------------------------------------------------------------

      Implicit    None

      Integer  :: n

      Real(l_) :: a(n)

      Integer  :: i

! ----------------------------------------------------------------------

      Do i = 1, n

         a(i) = dran0()

      End Do

! ----------------------------------------------------------------------

      End Subroutine ranfil

      Subroutine rinit

! ----------------------------------------------------------------------

      Implicit    None

      Integer  :: i

! ----------------------------------------------------------------------

      Do i = 1, 97                               ! --- Warming up phase.

         x1 = Mod( a1*x1 + c1, m1 )

         x2 = Mod( a2*x2 + c2, m2 )

      End Do

! ----------------------------------------------------------------------

      End Subroutine rinit

! -----------------------------------------------------------------

      End Module ran_module