MODULE randdp
! Replaces COMMON /randpp/ and defines dp

IMPLICIT NONE
INTEGER, PARAMETER, PUBLIC :: dp = SELECTED_REAL_KIND(15, 60)

REAL (dp), SAVE, PUBLIC    :: poly(101), other, offset
INTEGER, SAVE, PUBLIC      :: index

CONTAINS

!   Nick Maclaren's double precision random number generator.
!   This version, which is compatible with Lahey's ELF90 compiler,
!   is by Alan Miller ( alan @ vic.cmis.csiro.au, www.vic.cmis.csiro.au/~alan )

!   Latest revision - 18 December 1997

!   Copyright (C) 1992  N.M. Maclaren
!   Copyright (C) 1992  The University of Cambridge

!   This software may be reproduced and used freely, provided that all
!   users of it agree that the copyright holders are not liable for any
!   damage or injury caused by use of this software and that this
!   condition is passed onto all subsequent recipients of the software,
!   whether modified or not.



SUBROUTINE sdprnd (iseed)

INTEGER, INTENT(IN)  :: iseed

! Local variables
REAL (dp)            :: x
REAL (dp), PARAMETER :: xmod = 1000009711.0_dp, ymod = 33554432.0_dp
INTEGER              :: ix, iy, iz, i
LOGICAL, SAVE        :: inital = .TRUE.

!   ISEED should be set to an integer between 0 and 9999 inclusive;
!   a value of 0 will initialise the generator only if it has not
!   already been done.

IF (inital .OR. iseed /= 0) THEN
  inital = .false.
ELSE
  RETURN
END IF

!   INDEX must be initialised to an integer between 1 and 101 inclusive,
!   POLY(1...N) to integers between 0 and 1000009710 inclusive (not all 0),
!   and OTHER to a non-negative proper fraction with denominator 33554432.
!   It uses the Wichmann-Hill generator to do this.

ix = MOD(ABS(iseed), 10000) + 1
iy = 2*ix + 1
iz = 3*ix + 1
DO i = -10,101
  IF (i >= 1) poly(i) = AINT(xmod*x)
  ix = MOD(171*ix, 30269)
  iy = MOD(172*iy, 30307)
  iz = MOD(170*iz, 30323)
  x = MOD(DBLE(ix)/30269.0_dp + DBLE(iy)/30307.0_dp + DBLE(iz)/30323.0_dp, 1.0_dp)
END DO
other = AINT(ymod*x)/ymod
offset = 1.0_dp/ymod
INDEX = 1

RETURN
END SUBROUTINE sdprnd



FUNCTION dprand() RESULT(fn_val)

REAL (dp)            :: fn_val

! Local variables
REAL (dp)            :: x, y

! N.B. ymod has been removed from the previous DATA statement; it caused a
!      fatal error as it is not used.
REAL (dp), PARAMETER :: xmod = 1000009711.0_dp, xmod2 = 2000019422.0_dp, &
                        xmod4 = 4000038844.0_dp, tiny = 1.0E-17_dp,  &
                        zero = 0.0_dp, one = 1.0_dp
INTEGER              :: n
LOGICAL, SAVE        :: inital = .TRUE.

!   This returns a uniform (0,1) random number, with extremely good
!   uniformity properties.  It assumes that REAL (dp) provides
!   at least 33 bits of accuracy, and uses a power of two base.

IF (inital) THEN
  CALL sdprnd (0)
  inital = .false.
END IF

!   See [Knuth] for why this implements the algorithm described in the paper.
!   Note that this code is tuned for machines with fast REAL (dp), but
!   slow multiply and divide; many, many other options are possible.

n = INDEX - 64
IF (n <= 0) n = n + 101
x = poly(INDEX) + poly(INDEX)
x = xmod4 - poly(n) - poly(n) - x - x - poly(INDEX)
IF (x < zero) THEN
  IF (x < -xmod) x = x + xmod2
  IF (x < zero) x = x + xmod
ELSE
  IF (x >= xmod2) THEN
    x = x - xmod2
    IF (x >= xmod) x = x - xmod
  END IF
  IF (x >= xmod) x = x - xmod
END IF
poly(INDEX) = x
INDEX = INDEX + 1
IF (INDEX > 101) INDEX = INDEX - 101

!   Add in the second generator modulo 1, and force to be non-zero.
!   The restricted ranges largely cancel themselves out.

DO
  y = 37.0_dp*other + offset
  other = y - AINT(y)
  IF (other /= zero) EXIT
END DO

x = x/xmod + other
IF (x >= one) x = x - one
fn_val = x + tiny

RETURN
END FUNCTION dprand

END MODULE randdp