      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
