PRACE NPT

      Function ran1(idum)                             Result( ran_out )

      Use      numerics

      Implicit None

      Integer  idum

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

! --- ran1 returns a uniform deviate in (0,1).

! --- The algorithm is taken from Press & Teukolsky et.al. and

!     based on the linear congruential method  with choices for

!     M, IA, and IC that are given by D. Knuth in "Semi-numerical

!     algorithms.

! --- Input-parameters:

!     Integer - idum. When idum < 0 the sequence of random values

!                     When idum >= 0, ran1 returns the next value

!                     in the sequence. When ran1 is called for

!                     the first time it is also initialised.

! --- Output-parameters:

!     Integer  - idum   : Next value of seed produced by ran1.

!     Real(l_) - ran_out: the random number returned.

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

      Real(l_)            :: r(97), ran_out

      Integer             :: iff, x1, x2, x3, j

!

! --- Definitions of the three linear congruences used in generating

!     the random number.

!

      Integer, Parameter  :: m1 = 259200, m2 = 134456, m3 = 243000,

     &                       a1 =   7141, a2 =   8121, a3 =   4561,

     &                       c1 =  54773, c2 =  28411, c3 =  51349

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

!

      Save                   iff, r, x1, x2, x3

      Data                   iff/0/

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

! --- (Re)initialise if required.

      If ( idum < 0 .OR. iff == 0 ) Then

         iff = 1

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

! --- Seed first generator.

         x1 = Mod( c1 - idum,  m1 )

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

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

! --- Use it to seed the second generator.

         x2 = Mod( x1, m2 )

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

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

! --- Use generator 1 again to seed generator 3.

         x3 = Mod( a1*x1, m3 )

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

! --- Now fill array with random values, using gen. 2 for the high

!     order bits and gen. 1 for the low order bits.

         Do j = 1,97

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

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

            r(j) = ( Real( x1, Kind = l_ ) +

     &               Real( x2, Kind = l_ ) * rm2 ) * rm1

         End Do

         idum = 1

      End If

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

! --- This section is only reached when no (re)initialisation takes

!     place. A new random number is generated to fill the place of

!     a randomly picked element from array 'r' (the selection of the

!     index is done by gen. 3).

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

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

      x3 = Mod( a3*x3 + c3, m3 )

      j = 1 + ( 97 + x3 )/m3

      ran_out = r(j)

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

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

      End Function ran1