BROWSER/out_doc80/mode__random_8F90_source.html

 !SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier

 !SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence

 !SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt

 !SFX_LIC for details. version 1.

 MODULE mode_random

 ! A module for random number generation from the following distributions:

 !

 !     Distribution        Function/subroutine name

 !

 !     Normal (Gaussian)   random_normal

 !     Gamma   random_gamma

 !     Chi-squared         random_chisq

 !     Exponential         random_exponential

 !     Weibull random_Weibull

 !     Beta    random_beta

 !     t       random_t

 !     Multivariate normal random_mvnorm

 !     Generalized inverse Gaussian    random_inv_gauss

 !     Poisson random_Poisson

 !     Binomialrandom_binomial1   *

 ! random_binomial2   *

 !     Negative binomial   random_neg_binomial

 !     von Mises           random_von_Mises

 !     Cauchy  random_Cauchy

 !

 !  Generate a random ordering of the integers 1 .. N

 ! random_order

 !     Initialize (seed) the uniform random number generator for ANY compiler

 ! seed_random_number


 !     Lognormal - see note below.


 !  ** Two functions are provided for the binomial distribution.

 !  If the parameter values remain constant, it is recommended that the

 !  first function is used (random_binomial1).   If one or both of the

 !  parameters change, use the second function (random_binomial2).


 ! The compilers own random number generator, SUBROUTINE RANDOM_NUMBER(r),

 ! is used to provide a source of uniformly distributed random numbers.


 ! N.B. At this stage, only one random number is generated at each call to

 !      one of the functions above.


 ! The module uses the following functions which are included here:

 ! bin_prob to calculate a single binomial probability

 ! lngamma  to calculate the logarithm to base e of the gamma function


 ! Some of the code is adapted from Dagpunar's book:

 !     Dagpunar, J. 'Principles of random variate generation'

 !     Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9

 !

 ! In most of Dagpunar's routines, there is a test to see whether the value

 ! of one or two floating-point parameters has changed since the last call.

 ! These tests have been replaced by using a logical variable FIRST.

 ! This should be set to .TRUE. on the first call using new values of the

 ! parameters, and .FALSE. if the parameter values are the same as for the

 ! previous call.


 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

 ! Lognormal distribution

 ! If X has a lognormal distribution, then log(X) is normally distributed.

 ! Here the logarithm is the natural logarithm, that is to base e, sometimes

 ! denoted as ln.  To generate random variates from this distribution, generate

 ! a random deviate from the normal distribution with mean and variance equal

 ! to the mean and variance of the logarithms of X, then take its exponential.


 ! Relationship between the mean & variance of log(X) and the mean & variance

 ! of X, when X has a lognormal distribution.

 ! Let m = mean of log(X), and s^2 = variance of log(X)

 ! Then

 ! mean of X     = exp(m + 0.5s^2)

 ! variance of X = (mean(X))^2.[exp(s^2) - 1]


 ! In the reverse direction (rarely used)

 ! variance of log(X) = log[1 + var(X)/(mean(X))^2]

 ! mean of log(X)     = log(mean(X) - 0.5var(log(X))


 ! N.B. The above formulae relate to population parameters; they will only be

 !      approximate if applied to sample values.

 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


 ! Version 1.13, 2 October 2000

 ! Changes from version 1.01

 ! 1. The random_order, random_Poisson & random_binomial routines have been

 !    replaced with more efficient routines.

 ! 2. A routine, seed_random_number, has been added to seed the uniform random

 !    number generator.   This requires input of the required number of seeds

 !    for the particular compiler from a specified I/O unit such as a keyboard.

 ! 3. Made compatible with Lahey's ELF90.

 ! 4. Marsaglia & Tsang algorithm used for random_gamma when shape parameter > 1.

 ! 5. INTENT for array f corrected in random_mvnorm.


 !     Author: Alan Miller

 ! CSIRO Division of Mathematical & Information Sciences

 ! Private Bag 10, Clayton South MDC

 ! Clayton 3169, Victoria, Australia

 !     Phone: (+61) 3 9545-8016      Fax: (+61) 3 9545-8080

 !     e-mail: amiller @ bigpond.net.au

 !

 USE yomhook   ,ONLY : lhook,   dr_hook

 USE parkind1  ,ONLY : jprb

 !

 IMPLICIT NONE

 !

 REAL :: XHALF = 0.5

 !REAL :: XZERO = 0.0, XONE = 1.0, XTWO = 2.0,   &

 !        XVSMALL = TINY(1.0), XVLARGE = HUGE(1.0)

 !INTEGER, PARAMETER :: IDP = SELECTED_REAL_KIND(12, 60)

 !

  CONTAINS

 !

 !

 SUBROUTINE init_random_seed()

 !

 INTEGER, DIMENSION(:), ALLOCATABLE :: iseed

 INTEGER :: i, iin, iclock

 REAL(KIND=JPRB) :: zhook_handle

 !

 IF (lhook) CALL dr_hook('MODE_RANDOM_NUMBER:INIT_RANDOM_SEED',0,zhook_handle)

 !

  CALL random_seed(size = iin)

 ALLOCATE(iseed(iin))

 !

  CALL system_clock(count=iclock)

 !

 iseed(:) = iclock + 37 * (/ (i - 1, i = 1, iin) /)

 !

  CALL random_seed(put = iseed)

 !

 DEALLOCATE(iseed)

 !

 IF (lhook) CALL dr_hook('MODE_RANDOM_NUMBER:INIT_RANDOM_SEED',1,zhook_handle)

 !

 END SUBROUTINE init_random_seed

 !

 !

 FUNCTION random_normal() RESULT(PFN_VAL)

 !

 ! Adapted from the following Fortran 77 code

 !      ALGORITHM 712, COLLECTED ALGORITHMS FROM ACM.

 !      THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,

 !      VOL. 18, NO. 4, DECEMBER, 1992, PP. 434-435.

 !

 !  The function random_normal() returns a normally distributed pseudo-random

 !  number with zero mean and unit variance.

 !

 !  The algorithm uses the ratio of uniforms method of A.J. Kinderman

 !  and J.F. Monahan augmented with quadratic bounding curves.

 !

 REAL :: pfn_val

 !

 !     Local variables

 REAL :: zs = 0.449871, zt = -0.386595, za = 0.19600, zb = 0.25472, &

         zr1 = 0.27597, zr2 = 0.27846, zu, zv, zx, zy, zq

 !

 REAL(KIND=JPRB) :: zhook_handle

 !

 !     Generate P = (u,v) uniform in rectangle enclosing acceptance region

 !

 IF (lhook) CALL dr_hook('MODE_RANDOM_NUMBER:RANDOM_NORMAL',0,zhook_handle)

 !

 DO

   !

   CALL random_number(zu)

   CALL random_number(zv)

   !

   zv = 1.7156 * (zv - xhalf)

   !

   !     Evaluate the quadratic form

   !

   zx = zu - zs

   zy = abs(zv) - zt

   zq = zx**2 + zy*(za*zy - zb*zx)

   !

   !     Accept P if inside inner ellipse

   IF (zq < zr1) EXIT

   !

   !     Reject P if outside outer ellipse

   IF (zq > zr2) cycle

   !

   !     Reject P if outside acceptance region

   IF (zv**2 < -4.0*log(zu)*zu**2) EXIT

   !

 END DO

 !

 !     Return ratio of P's coordinates as the normal deviate

 pfn_val = zv/zu

 !

 IF (lhook) CALL dr_hook('MODE_RANDOM_NUMBER:RANDOM_NORMAL',1,zhook_handle)

 !

 END FUNCTION random_normal

 !

 !FUNCTION random_gamma(s, first) RESULT(fn_val)

 !

 !! Adapted from Fortran 77 code from the book:

 !!     Dagpunar, J. 'Principles of random variate generation'

 !!     Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9

 !

 !!     FUNCTION GENERATES A RANDOM GAMMA VARIATE.

 !!     CALLS EITHER random_gamma1 (S > 1.0)

 !!     OR random_exponential (S = 1.0)

 !!     OR random_gamma2 (S < 1.0).

 !

 !!     S = SHAPE PARAMETER OF DISTRIBUTION (0 < REAL).

 !

 !REAL, INTENT(IN)    :: s

 !LOGICAL, INTENT(IN) :: first

 !REAL    :: fn_val

 !

 !IF (s <= XZERO) THEN

 !  WRITE(*, *) 'SHAPE PARAMETER VALUE MUST BE POSITIVE'

 !  STOP

 !END IF

 !

 !IF (s > XONE) THEN

 !  fn_val = random_gamma1(s, first)

 !ELSE IF (s < XONE) THEN

 !  fn_val = random_gamma2(s, first)

 !ELSE

 !  fn_val = random_expXONEntial()

 !END IF

 !

 !RETURN

 !END FUNCTION random_gamma

 !

 !

 !

 !FUNCTION random_gamma1(s, first) RESULT(fn_val)

 !

 !! Uses the algorithm in

 !! Marsaglia, G. and Tsang, W.W. (2000) `A simple method for generating

 !! gamma variables', Trans. om Math. Software (TOMS), vol.26(3), pp.363-372.

 !

 !! Generates a random gamma deviate for shape parameter s >= 1.

 !

 !REAL, INTENT(IN)    :: s

 !LOGICAL, INTENT(IN) :: first

 !REAL    :: fn_val

 !

 !! Local variables

 !REAL, SAVE  :: c, d

 !REAL        :: u, v, x

 !

 !IF (first) THEN

 !  d = s - XONE/3.

 !  c = XONE/SQRT(9.0*d)

 !END IF

 !

 !! Start of main loop

 !DO

 !

 !! Generate v = (1+cx)^3 where x is random normal; repeat if v <= 0.

 !

 !  DO

 !    x = random_normal()

 !    v = (XONE + c*x)**3

 !    IF (v > XZERO) EXIT

 !  END DO

 !

 !! Generate uniform variable U

 !

 !  CALL RANDOM_NUMBER(u)

 !  IF (u < XONE - 0.0331*x**4) THEN

 !    fn_val = d*v

 !    EXIT

 !  ELSE IF (LOG(u) < XHALF*x**2 + d*(XONE - v + LOG(v))) THEN

 !    fn_val = d*v

 !    EXIT

 !  END IF

 !END DO

 !

 !RETURN

 !END FUNCTION random_gamma1

 !

 !

 !

 !FUNCTION random_gamma2(s, first) RESULT(fn_val)

 !

 !! Adapted from Fortran 77 code from the book:

 !!     Dagpunar, J. 'Principles of random variate generation'

 !!     Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9

 !

 !! FUNCTION GENERATES A RANDOM VARIATE IN [0,INFINITY) FROM

 !! A GAMMA DISTRIBUTION WITH DENSITY PROPORTIONAL TO

 !! GAMMA2**(S-1) * EXP(-GAMMA2),

 !! USING A SWITCHING METHOD.

 !

 !!    S = SHAPE PARAMETER OF DISTRIBUTION

 !!          (REAL < 1.0)

 !

 !REAL, INTENT(IN)    :: s

 !LOGICAL, INTENT(IN) :: first

 !REAL    :: fn_val

 !

 !!     Local variables

 !REAL       :: r, x, w

 !REAL, SAVE :: a, p, c, uf, vr, d

 !

 !IF (s <= XZERO .OR. s >= XONE) THEN

 !  WRITE(*, *) 'SHAPE PARAMETER VALUE OUTSIDE PERMITTED RANGE'

 !  STOP

 !END IF

 !

 !IF (first) THEN! Initialization, if necessary

 !  a = XONE - s

 !  p = a/(a + s*EXP(-a))

 !  IF (s < XVSMALL) THEN

 !    WRITE(*, *) 'SHAPE PARAMETER VALUE TOO SMALL'

 !    STOP

 !  END IF

 !  c = XONE/s

 !  uf = p*(XVSMALL/a)**s

 !  vr = XONE - XVSMALL

 !  d = a*LOG(a)

 !END IF

 !

 !DO

 !  CALL RANDOM_NUMBER(r)

 !  IF (r >= vr) THEN

 !    CYCLE

 !  ELSE IF (r > p) THEN

 !    x = a - LOG((XONE - r)/(XONE - p))

 !    w = a*LOG(x)-d

 !  ELSE IF (r > uf) THEN

 !    x = a*(r/p)**c

 !    w = x

 !  ELSE

 !    fn_val = XZERO

 !    RETURN

 !  END IF

 !

 !  CALL RANDOM_NUMBER(r)

 !  IF (XONE-r <= w .AND. r > XZERO) THEN

 !    IF (r*(w + XONE) >= XONE) CYCLE

 !    IF (-LOG(r) <= w) CYCLE

 !  END IF

 !  EXIT

 !END DO

 !

 !fn_val = x

 !RETURN

 !

 !END FUNCTION random_gamma2

 !

 !

 !

 !FUNCTION random_chisq(ndf, first) RESULT(fn_val)

 !

 !!     Generates a random variate from the chi-squared distribution with

 !!     ndf degrees of freedom

 !

 !INTEGER, INTENT(IN) :: ndf

 !LOGICAL, INTENT(IN) :: first

 !REAL    :: fn_val

 !

 !fn_val = XTWO * random_gamma(XHALF*ndf, first)

 !RETURN

 !

 !END FUNCTION random_chisq

 !

 !!

 !FUNCTION random_exponential() RESULT(fn_val)

 !

 !! Adapted from Fortran 77 code from the book:

 !!     Dagpunar, J. 'Principles of random variate generation'

 !!     Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9

 !

 !! FUNCTION GENERATES A RANDOM VARIATE IN [0,INFINITY) FROM

 !! A NEGATIVE EXPONENTIAL DlSTRIBUTION WlTH DENSITY PROPORTIONAL

 !! TO EXP(-random_exponential), USING INVERSION.

 !

 !REAL  :: fn_val

 !

 !!     Local variable

 !REAL  :: r

 !

 !DO

 !  CALL RANDOM_NUMBER(r)

 !  IF (r > XZERO) EXIT

 !END DO

 !

 !fn_val = -LOG(r)

 !RETURN

 !

 !END FUNCTION random_exponential

 !

 !

 !

 !FUNCTION random_Weibull(a) RESULT(fn_val)

 !

 !!     Generates a random variate from the Weibull distribution with

 !!     probability density:

 !!          a

 !!   a-1  -x

 !!     f(x) = a.x    e

 !

 !REAL, INTENT(IN) :: a

 !REAL :: fn_val

 !

 !!     For speed, there is no checking that a is not zero or very small.

 !

 !fn_val = random_exponential() ** (XONE/a)

 !RETURN

 !

 !END FUNCTION random_Weibull

 !

 !

 !

 !FUNCTION random_beta(aa, bb, first) RESULT(fn_val)

 !

 !! Adapted from Fortran 77 code from the book:

 !!     Dagpunar, J. 'Principles of random variate generation'

 !!     Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9

 !

 !! FUNCTION GENERATES A RANDOM VARIATE IN [0,1]

 !! FROM A BETA DISTRIBUTION WITH DENSITY

 !! PROPORTIONAL TO BETA**(AA-1) * (1-BETA)**(BB-1).

 !! USING CHENG'S LOG LOGISTIC METHOD.

 !

 !!     AA = SHAPE PARAMETER FROM DISTRIBUTION (0 < REAL)

 !!     BB = SHAPE PARAMETER FROM DISTRIBUTION (0 < REAL)

 !

 !REAL, INTENT(IN)    :: aa, bb

 !LOGICAL, INTENT(IN) :: first

 !REAL    :: fn_val

 !

 !!     Local variables

 !REAL, PARAMETER  :: aln4 = 1.3862944

 !REAL :: a, b, g, r, s, x, y, z

 !REAL, SAVE       :: d, f, h, t, c

 !LOGICAL, SAVE    :: swap

 !

 !IF (aa <= XZERO .OR. bb <= XZERO) THEN

 !  WRITE(*, *) 'IMPERMISSIBLE SHAPE PARAMETER VALUE(S)'

 !  STOP

 !END IF

 !

 !IF (first) THEN! Initialization, if necessary

 !  a = aa

 !  b = bb

 !  swap = b > a

 !  IF (swap) THEN

 !    g = b

 !    b = a

 !    a = g

 !  END IF

 !  d = a/b

 !  f = a+b

 !  IF (b > XONE) THEN

 !    h = SQRT((XTWO*a*b - f)/(f - XTWO))

 !    t = XONE

 !  ELSE

 !    h = b

 !    t = XONE/(XONE + (a/(XVLARGE*b))**b)

 !  END IF

 !  c = a+h

 !END IF

 !

 !DO

 !  CALL RANDOM_NUMBER(r)

 !  CALL RANDOM_NUMBER(x)

 !  s = r*r*x

 !  IF (r < XVSMALL .OR. s <= XZERO) CYCLE

 !  IF (r < t) THEN

 !    x = LOG(r/(XONE - r))/h

 !    y = d*EXP(x)

 !    z = c*x + f*LOG((XONE + d)/(XONE + y)) - aln4

 !    IF (s - XONE > z) THEN

 !      IF (s - s*z > XONE) CYCLE

 !      IF (LOG(s) > z) CYCLE

 !    END IF

 !    fn_val = y/(XONE + y)

 !  ELSE

 !    IF (4.0*s > (XONE + XONE/d)**f) CYCLE

 !    fn_val = XONE

 !  END IF

 !  EXIT

 !END DO

 !

 !IF (swap) fn_val = XONE - fn_val

 !RETURN

 !END FUNCTION random_beta

 !

 !

 !

 !FUNCTION random_t(m) RESULT(fn_val)

 !

 !! Adapted from Fortran 77 code from the book:

 !!     Dagpunar, J. 'Principles of random variate generation'

 !!     Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9

 !

 !! FUNCTION GENERATES A RANDOM VARIATE FROM A

 !! T DISTRIBUTION USING KINDERMAN AND MONAHAN'S RATIO METHOD.

 !

 !!     M = DEGREES OF FREEDOM OF DISTRIBUTION

 !!           (1 <= 1NTEGER)

 !

 !INTEGER, INTENT(IN) :: m

 !REAL    :: fn_val

 !

 !!     Local variables

 !REAL, SAVE      :: s, c, a, f, g

 !REAL:: r, x, v

 !

 !REAL, PARAMETER :: three = 3.0, four = 4.0, quart = 0.25,   &

 !       five = 5.0, sixteen = 16.0

 !INTEGER         :: mm = 0

 !

 !IF (m < 1) THEN

 !  WRITE(*, *) 'IMPERMISSIBLE DEGREES OF FREEDOM'

 !  STOP

 !END IF

 !

 !IF (m /= mm) THEN        ! Initialization, if necessary

 !  s = m

 !  c = -quart*(s + XONE)

 !  a = four/(XONE + XONE/s)**c

 !  f = sixteen/a

 !  IF (m > 1) THEN

 !    g = s - XONE

 !    g = ((s + XONE)/g)**c*SQRT((s+s)/g)

 !  ELSE

 !    g = XONE

 !  END IF

 !  mm = m

 !END IF

 !

 !DO

 !  CALL RANDOM_NUMBER(r)

 !  IF (r <= XZERO) CYCLE

 !  CALL RANDOM_NUMBER(v)

 !  x = (XTWO*v - XONE)*g/r

 !  v = x*x

 !  IF (v > five - a*r) THEN

 !    IF (m >= 1 .AND. r*(v + three) > f) CYCLE

 !    IF (r > (XONE + v/s)**c) CYCLE

 !  END IF

 !  EXIT

 !END DO

 !

 !fn_val = x

 !RETURN

 !END FUNCTION random_t

 !

 !

 !

 !SUBROUTINE random_mvnorm(n, h, d, f, first, x, ier)

 !

 !! Adapted from Fortran 77 code from the book:

 !!     Dagpunar, J. 'Principles of random variate generation'

 !!     Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9

 !

 !! N.B. An extra argument, ier, has been added to Dagpunar's routine

 !

 !!     SUBROUTINE GENERATES AN N VARIATE RANDOM NORMAL

 !!     VECTOR USING A CHOLESKY DECOMPOSITION.

 !

 !! ARGUMENTS:

 !!        N = NUMBER OF VARIATES IN VECTOR

 !!           (INPUT,INTEGER >= 1)

 !!     H(J) = J'TH ELEMENT OF VECTOR OF MEANS

 !!           (INPUT,REAL)

 !!     X(J) = J'TH ELEMENT OF DELIVERED VECTOR

 !!           (OUTPUT,REAL)

 !!

 !!    D(J*(J-1)/2+I) = (I,J)'TH ELEMENT OF VARIANCE MATRIX (J> = I)

 !!(INPUT,REAL)

 !!    F((J-1)*(2*N-J)/2+I) = (I,J)'TH ELEMENT OF LOWER TRIANGULAR

 !!           DECOMPOSITION OF VARIANCE MATRIX (J <= I)

 !!(OUTPUT,REAL)

 !

 !!    FIRST = .TRUE. IF THIS IS THE FIRST CALL OF THE ROUTINE

 !!    OR IF THE DISTRIBUTION HAS CHANGED SINCE THE LAST CALL OF THE ROUTINE.

 !!    OTHERWISE SET TO .FALSE.

 !!(INPUT,LOGICAL)

 !

 !!    ier = 1 if the input covariance matrix is not +ve definite

 !!        = 0 otherwise

 !

 !INTEGER, INTENT(IN)   :: n

 !REAL, INTENT(IN)      :: h(:), d(:)   ! d(n*(n+1)/2)

 !REAL, INTENT(IN OUT)  :: f(:)         ! f(n*(n+1)/2)

 !REAL, INTENT(OUT)     :: x(:)

 !LOGICAL, INTENT(IN)   :: first

 !INTEGER, INTENT(OUT)  :: ier

 !

 !!     Local variables

 !INTEGER       :: j, i, m

 !REAL          :: y, v

 !INTEGER, SAVE :: n2

 !

 !IF (n < 1) THEN

 !  WRITE(*, *) 'SIZE OF VECTOR IS NON POSITIVE'

 !  STOP

 !END IF

 !

 !ier = 0

 !IF (first) THEN! Initialization, if necessary

 !  n2 = 2*n

 !  IF (d(1) < XZERO) THEN

 !    ier = 1

 !    RETURN

 !  END IF

 !

 !  f(1) = SQRT(d(1))

 !  y = XONE/f(1)

 !  DO j = 2,n

 !    f(j) = d(1+j*(j-1)/2) * y

 !  END DO

 !

 !  DO i = 2,n

 !    v = d(i*(i-1)/2+i)

 !    DO m = 1,i-1

 !      v = v - f((m-1)*(n2-m)/2+i)**2

 !    END DO

 !

 !    IF (v < XZERO) THEN

 !      ier = 1

 !      RETURN

 !    END IF

 !

 !    v = SQRT(v)

 !    y = XONE/v

 !    f((i-1)*(n2-i)/2+i) = v

 !    DO j = i+1,n

 !      v = d(j*(j-1)/2+i)

 !      DO m = 1,i-1

 !        v = v - f((m-1)*(n2-m)/2+i)*f((m-1)*(n2-m)/2 + j)

 !      END DO ! m = 1,i-1

 !      f((i-1)*(n2-i)/2 + j) = v*y

 !    END DO ! j = i+1,n

 !  END DO ! i = 2,n

 !END IF

 !

 !x(1:n) = h(1:n)

 !DO j = 1,n

 !  y = random_normal()

 !  DO i = j,n

 !    x(i) = x(i) + f((j-1)*(n2-j)/2 + i) * y

 !  END DO ! i = j,n

 !END DO ! j = 1,n

 !

 !RETURN

 !END SUBROUTINE random_mvnorm

 !

 !

 !

 !FUNCTION random_inv_gauss(h, b, first) RESULT(fn_val)

 !

 !! Adapted from Fortran 77 code from the book:

 !!     Dagpunar, J. 'Principles of random variate generation'

 !!     Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9

 !

 !! FUNCTION GENERATES A RANDOM VARIATE IN [0,INFINITY] FROM

 !! A REPARAMETERISED GENERALISED INVERSE GAUSSIAN (GIG) DISTRIBUTION

 !! WITH DENSITY PROPORTIONAL TO  GIG**(H-1) * EXP(-0.5*B*(GIG+1/GIG))

 !! USING A RATIO METHOD.

 !

 !!     H = PARAMETER OF DISTRIBUTION (0 <= REAL)

 !!     B = PARAMETER OF DISTRIBUTION (0 < REAL)

 !

 !REAL, INTENT(IN)    :: h, b

 !LOGICAL, INTENT(IN) :: first

 !REAL    :: fn_val

 !

 !!     Local variables

 !REAL:: ym, xm, r, w, r1, r2, x

 !REAL, SAVE      :: a, c, d, e

 !REAL, PARAMETER :: quart = 0.25

 !

 !IF (h < XZERO .OR. b <= XZERO) THEN

 !  WRITE(*, *) 'IMPERMISSIBLE DISTRIBUTION PARAMETER VALUES'

 !  STOP

 !END IF

 !

 !IF (first) THEN! Initialization, if necessary

 !  IF (h > quart*b*SQRT(XVLARGE)) THEN

 !    WRITE(*, *) 'THE RATIO H:B IS TOO SMALL'

 !    STOP

 !  END IF

 !  e = b*b

 !  d = h + XONE

 !  ym = (-d + SQRT(d*d + e))/b

 !  IF (ym < XVSMALL) THEN

 !    WRITE(*, *) 'THE VALUE OF B IS TOO SMALL'

 !    STOP

 !  END IF

 !

 !  d = h - XONE

 !  xm = (d + SQRT(d*d + e))/b

 !  d = XHALF*d

 !  e = -quart*b

 !  r = xm + XONE/xm

 !  w = xm*ym

 !  a = w**(-XHALF*h) * SQRT(xm/ym) * EXP(-e*(r - ym - XONE/ym))

 !  IF (a < XVSMALL) THEN

 !    WRITE(*, *) 'THE VALUE OF H IS TOO LARGE'

 !    STOP

 !  END IF

 !  c = -d*LOG(xm) - e*r

 !END IF

 !

 !DO

 !  CALL RANDOM_NUMBER(r1)

 !  IF (r1 <= XZERO) CYCLE

 !  CALL RANDOM_NUMBER(r2)

 !  x = a*r2/r1

 !  IF (x <= XZERO) CYCLE

 !  IF (LOG(r1) < d*LOG(x) + e*(x + XONE/x) + c) EXIT

 !END DO

 !

 !fn_val = x

 !

 !RETURN

 !END FUNCTION random_inv_gauss

 !

 !

 !

 !FUNCTION random_Poisson(mu, first) RESULT(ival)

 !!**********************************************************************

 !!     Translated to Fortran 90 by Alan Miller from:

 !!   RANLIB

 !!

 !!     Library of Fortran Routines for Random Number Generation

 !!

 !!        Compiled and Written by:

 !!

 !! Barry W. Brown

 !!  James Lovato

 !!

 !! Department of Biomathematics, Box 237

 !! The University of Texas, M.D. Anderson Cancer Center

 !! 1515 Holcombe Boulevard

 !! Houston, TX      77030

 !!

 !! This work was supported by grant CA-16672 from the National Cancer Institute.

 !

 !!        GENerate POIsson random deviate

 !

 !!    Function

 !

 !! Generates a single random deviate from a Poisson distribution with mean mu.

 !

 !!    Arguments

 !

 !!     mu --> The mean of the Poisson distribution from which

 !!a random deviate is to be generated.

 !!      REAL mu

 !

 !!      Method

 !

 !!     For details see:

 !

 !!   Ahrens, J.H. and Dieter, U.

 !!   Computer Generation of Poisson Deviates

 !!   From Modified Normal Distributions.

 !!   ACM Trans. Math. Software, 8, 2

 !!   (June 1982),163-179

 !

 !!     TABLES: COEFFICIENTS A0-A7 FOR STEP F. FACTORIALS FACT

 !!     COEFFICIENTS A(K) - FOR PX = FK*V*V*SUM(A(K)*V**K)-DEL

 !

 !!     SEPARATION OF CASES A AND B

 !

 !!     .. Scalar Arguments ..

 !REAL, INTENT(IN)    :: mu

 !LOGICAL, INTENT(IN) :: first

 !INTEGER :: ival

 !!     ..

 !!     .. Local Scalars ..

 !REAL          :: b1, b2, c, c0, c1, c2, c3, del, difmuk, e, fk, fx, fy, g,  &

 !     omega, px, py, t, u, v, x, xx

 !REAL, SAVE    :: s, d, p, q, p0

 !INTEGER       :: j, k, kflag

 !LOGICAL, SAVE :: full_init

 !INTEGER, SAVE :: l, m

 !!     ..

 !!     .. Local Arrays ..

 !REAL, SAVE    :: pp(35)

 !!     ..

 !!     .. Data statements ..

 !REAL, PARAMETER :: a0 = -.5, a1 = .3333333, a2 = -.2500068, a3 = .2000118,  &

 !       a4 = -.1661269, a5 = .1421878, a6 = -.1384794,   &

 !       a7 = .1250060

 !

 !REAL, PARAMETER :: fact(10) = (/ 1., 1., 2., 6., 24., 120., 720., 5040.,  &

 !         40320., 362880. /)

 !

 !!     ..

 !!     .. Executable Statements ..

 !IF (mu > 10.0) THEN

 !!     C A S E  A. (RECALCULATION OF S, D, L IF MU HAS CHANGED)

 !

 !  IF (first) THEN

 !    s = SQRT(mu)

 !    d = 6.0*mu*mu

 !

 !! THE POISSON PROBABILITIES PK EXCEED THE DISCRETE NORMAL

 !! PROBABILITIES FK WHENEVER K >= M(MU). L=IFIX(MU-1.1484)

 !! IS AN UPPER BOUND TO M(MU) FOR ALL MU >= 10 .

 !

 !    l = mu - 1.1484

 !    full_init = .false.

 !  END IF

 !

 !

 !!     STEP N. NORMAL SAMPLE - random_normal() FOR STANDARD NORMAL DEVIATE

 !

 !  g = mu + s*random_normal()

 !  IF (g > 0.0) THEN

 !    ival = g

 !

 !!     STEP I. IMMEDIATE ACCEPTANCE IF ival IS LARGE ENOUGH

 !

 !    IF (ival>=l) RETURN

 !

 !!     STEP S. SQUEEZE ACCEPTANCE - SAMPLE U

 !

 !    fk = ival

 !    difmuk = mu - fk

 !    CALL RANDOM_NUMBER(u)

 !    IF (d*u >= difmuk*difmuk*difmuk) RETURN

 !  END IF

 !

 !!     STEP P. PREPARATIONS FOR STEPS Q AND H.

 !! (RECALCULATIONS OF PARAMETERS IF NECESSARY)

 !! .3989423=(2*PI)**(-.5)  .416667E-1=1./24.  .1428571=1./7.

 !! THE QUANTITIES B1, B2, C3, C2, C1, C0 ARE FOR THE HERMITE

 !! APPROXIMATIONS TO THE DISCRETE NORMAL PROBABILITIES FK.

 !! C=.1069/MU GUARANTEES MAJORIZATION BY THE 'HAT'-FUNCTION.

 !

 !  IF (.NOT. full_init) THEN

 !    omega = .3989423/s

 !    b1 = .4166667E-1/mu

 !    b2 = .3*b1*b1

 !    c3 = .1428571*b1*b2

 !    c2 = b2 - 15.*c3

 !    c1 = b1 - 6.*b2 + 45.*c3

 !    c0 = 1. - b1 + 3.*b2 - 15.*c3

 !    c = .1069/mu

 !    full_init = .true.

 !  END IF

 !

 !  IF (g < 0.0) GO TO 50

 !

 !! 'SUBROUTINE' F IS CALLED (KFLAG=0 FOR CORRECT RETURN)

 !

 !  kflag = 0

 !  GO TO 70

 !

 !!     STEP Q. QUOTIENT ACCEPTANCE (RARE CASE)

 !

 !  40 IF (fy-u*fy <= py*EXP(px-fx)) RETURN

 !

 !!     STEP E. EXPONENTIAL SAMPLE - random_exponential() FOR STANDARD EXPONENTIAL

 !! DEVIATE E AND SAMPLE T FROM THE LAPLACE 'HAT'

 !! (IF T <= -.6744 THEN PK < FK FOR ALL MU >= 10.)

 !

 !  50 e = random_exponential()

 !  CALL RANDOM_NUMBER(u)

 !  u = u + u - XONE

 !  t = 1.8 + SIGN(e, u)

 !  IF (t <= (-.6744)) GO TO 50

 !  ival = mu + s*t

 !  fk = ival

 !  difmuk = mu - fk

 !

 !! 'SUBROUTINE' F IS CALLED (KFLAG=1 FOR CORRECT RETURN)

 !

 !  kflag = 1

 !  GO TO 70

 !

 !!     STEP H. HAT ACCEPTANCE (E IS REPEATED ON REJECTION)

 !

 !  60 IF (c*ABS(u) > py*EXP(px+e) - fy*EXP(fx+e)) GO TO 50

 !  RETURN

 !

 !!     STEP F. 'SUBROUTINE' F. CALCULATION OF PX, PY, FX, FY.

 !! CASE ival < 10 USES FACTORIALS FROM TABLE FACT

 !

 !  70 IF (ival>=10) GO TO 80

 !  px = -mu

 !  py = mu**ival/fact(ival+1)

 !  GO TO 110

 !

 !! CASE ival >= 10 USES POLYNOMIAL APPROXIMATION

 !! A0-A7 FOR ACCURACY WHEN ADVISABLE

 !! .8333333E-1=1./12.  .3989423=(2*PI)**(-.5)

 !

 !  80 del = .8333333E-1/fk

 !  del = del - 4.8*del*del*del

 !  v = difmuk/fk

 !  IF (ABS(v)>0.25) THEN

 !    px = fk*LOG(XONE + v) - difmuk - del

 !  ELSE

 !    px = fk*v*v* (((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v+a0) - del

 !  END IF

 !  py = .3989423/SQRT(fk)

 !  110 x = (XHALF - difmuk)/s

 !  xx = x*x

 !  fx = -XHALF*xx

 !  fy = omega* (((c3*xx + c2)*xx + c1)*xx + c0)

 !  IF (kflag <= 0) GO TO 40

 !  GO TO 60

 !

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

 !!     C A S E  B.    mu < 10

 !!     START NEW TABLE AND CALCULATE P0 IF NECESSARY

 !

 !ELSE

 !  IF (first) THEN

 !    m = MAX(1, INT(mu))

 !    l = 0

 !    p = EXP(-mu)

 !    q = p

 !    p0 = p

 !  END IF

 !

 !!     STEP U. UNIFORM SAMPLE FOR INVERSION METHOD

 !

 !  DO

 !    CALL RANDOM_NUMBER(u)

 !    ival = 0

 !    IF (u <= p0) RETURN

 !

 !!     STEP T. TABLE COMPARISON UNTIL THE END PP(L) OF THE

 !! PP-TABLE OF CUMULATIVE POISSON PROBABILITIES

 !! (0.458=PP(9) FOR MU=10)

 !

 !    IF (l == 0) GO TO 150

 !    j = 1

 !    IF (u > 0.458) j = MIN(l, m)

 !    DO k = j, l

 !      IF (u <= pp(k)) GO TO 180

 !    END DO

 !    IF (l == 35) CYCLE

 !

 !!     STEP C. CREATION OF NEW POISSON PROBABILITIES P

 !! AND THEIR CUMULATIVES Q=PP(K)

 !

 !    150 l = l + 1

 !    DO k = l, 35

 !      p = p*mu / k

 !      q = q + p

 !      pp(k) = q

 !      IF (u <= q) GO TO 170

 !    END DO

 !    l = 35

 !  END DO

 !

 !  170 l = k

 !  180 ival = k

 !  RETURN

 !END IF

 !

 !RETURN

 !END FUNCTION random_Poisson

 !

 !

 !

 !FUNCTION random_binomial1(n, p, first) RESULT(ival)

 !

 !! FUNCTION GENERATES A RANDOM BINOMIAL VARIATE USING C.D.Kemp's method.

 !! This algorithm is suitable when many random variates are required

 !! with the SAME parameter values for n & p.

 !

 !!    P = BERNOULLI SUCCESS PROBABILITY

 !!           (0 <= REAL <= 1)

 !!    N = NUMBER OF BERNOULLI TRIALS

 !!           (1 <= INTEGER)

 !!    FIRST = .TRUE. for the first call using the current parameter values

 !!          = .FALSE. if the values of (n,p) are unchanged from last call

 !

 !! Reference: Kemp, C.D. (1986). `A modal method for generating binomial

 !!variables', Commun. Statist. - Theor. Meth. 15(3), 805-813.

 !

 !INTEGER, INTENT(IN) :: n

 !REAL, INTENT(IN)    :: p

 !LOGICAL, INTENT(IN) :: first

 !INTEGER :: ival

 !

 !!     Local variables

 !

 !INTEGER         :: ru, rd

 !INTEGER, SAVE   :: r0

 !REAL:: u, pd, pu

 !REAL, SAVE      :: odds_ratio, p_r

 !REAL, PARAMETER :: ZZERO = 0.0, ZONE = 1.0

 !

 !IF (first) THEN

 !  r0 = (n+1)*p

 !  p_r = bin_prob(n, p, r0)

 !  odds_ratio = p / (ZONE - p)

 !END IF

 !

 !CALL RANDOM_NUMBER(u)

 !u = u - p_r

 !IF (u < ZZERO) THEN

 !  ival = r0

 !  RETURN

 !END IF

 !

 !pu = p_r

 !ru = r0

 !pd = p_r

 !rd = r0

 !DO

 !  rd = rd - 1

 !  IF (rd >= 0) THEN

 !    pd = pd * (rd+1) / (odds_ratio * (n-rd))

 !    u = u - pd

 !    IF (u < ZZERO) THEN

 !      ival = rd

 !      RETURN

 !    END IF

 !  END IF

 !

 !  ru = ru + 1

 !  IF (ru <= n) THEN

 !    pu = pu * (n-ru+1) * odds_ratio / ru

 !    u = u - pu

 !    IF (u < ZZERO) THEN

 !      ival = ru

 !      RETURN

 !    END IF

 !  END IF

 !END DO

 !

 !!     This point should not be reached, but just in case:

 !

 !ival = r0

 !RETURN

 !

 !END FUNCTION random_binomial1

 !

 !

 !

 !FUNCTION bin_prob(n, p, r) RESULT(fn_val)

 !!     Calculate a binomial probability

 !

 !INTEGER, INTENT(IN) :: n, r

 !REAL, INTENT(IN)    :: p

 !REAL    :: fn_val

 !

 !!     Local variable

 !REAL    :: ZONE = 1.0

 !

 !fn_val = EXP( lngamma(DBLE(n+1)) - lngamma(DBLE(r+1)) - lngamma(DBLE(n-r+1)) &

 !  + r*LOG(p) + (n-r)*LOG(ZONE - p) )

 !RETURN

 !

 !END FUNCTION bin_prob

 !

 !

 !

 !FUNCTION lngamma(x) RESULT(fn_val)

 !

 !! Logarithm to base e of the gamma function.

 !!

 !! Accurate to about 1.e-14.

 !! Programmer: Alan Miller

 !

 !! Latest revision of Fortran 77 version - 28 February 1988

 !

 !REAL (IDP), INTENT(IN) :: x

 !REAL (IDP) :: fn_val

 !

 !!       Local variables

 !

 !REAL (IDP) :: a1 = -4.166666666554424D-02, a2 = 2.430554511376954D-03,  &

 ! a3 = -7.685928044064347D-04, a4 = 5.660478426014386D-04,  &

 ! temp, arg, product, lnrt2pi = 9.189385332046727D-1,       &

 ! pi = 3.141592653589793D0

 !LOGICAL   :: reflect

 !

 !!       lngamma is not defined if x = 0 or a negative integer.

 !

 !IF (x > 0.d0) GO TO 10

 !IF (x /= INT(x)) GO TO 10

 !fn_val = 0.d0

 !RETURN

 !

 !!       If x < 0, use the reflection formula:

 !!   gamma(x) * gamma(1-x) = pi * cosec(pi.x)

 !

 !10 reflect = (x < 0.d0)

 !IF (reflect) THEN

 !  arg = 1.d0 - x

 !ELSE

 !  arg = x

 !END IF

 !

 !!       Increase the argument, if necessary, to make it > 10.

 !

 !product = 1.d0

 !20 IF (arg <= 10.d0) THEN

 !  product = product * arg

 !  arg = arg + 1.d0

 !  GO TO 20

 !END IF

 !

 !!  Use a polynomial approximation to Stirling's formula.

 !!  N.B. The real Stirling's formula is used here, not the simpler, but less

 !!       accurate formula given by De Moivre in a letter to Stirling, which

 !!       is the one usually quoted.

 !

 !arg = arg - 0.5D0

 !temp = 1.d0/arg**2

 !fn_val = lnrt2pi + arg * (LOG(arg) - 1.d0 + &

 !      (((a4*temp + a3)*temp + a2)*temp + a1)*temp) - LOG(product)

 !IF (reflect) THEN

 !  temp = SIN(pi * x)

 !  fn_val = LOG(pi/temp) - fn_val

 !END IF

 !RETURN

 !END FUNCTION lngamma

 !

 !

 !

 !FUNCTION random_binomial2(n, pp, first) RESULT(ival)

 !!**********************************************************************

 !!     Translated to Fortran 90 by Alan Miller from:

 !!      RANLIB

 !!

 !!     Library of Fortran Routines for Random Number Generation

 !!

 !!          Compiled and Written by:

 !!

 !!   Barry W. Brown

 !!    James Lovato

 !!

 !!   Department of Biomathematics, Box 237

 !!   The University of Texas, M.D. Anderson Cancer Center

 !!   1515 Holcombe Boulevard

 !!   Houston, TX      77030

 !!

 !! This work was supported by grant CA-16672 from the National Cancer Institute.

 !

 !!        GENerate BINomial random deviate

 !

 !!      Function

 !

 !!     Generates a single random deviate from a binomial

 !!     distribution whose number of trials is N and whose

 !!     probability of an event in each trial is P.

 !

 !!      Arguments

 !

 !!     N  --> The number of trials in the binomial distribution

 !!from which a random deviate is to be generated.

 !!      INTEGER N

 !

 !!     P  --> The probability of an event in each trial of the

 !!binomial distribution from which a random deviate

 !!is to be generated.

 !!      REAL P

 !

 !!     FIRST --> Set FIRST = .TRUE. for the first call to perform initialization

 !!   the set FIRST = .FALSE. for further calls using the same pair

 !!   of parameter values (N, P).

 !!      LOGICAL FIRST

 !

 !!     random_binomial2 <-- A random deviate yielding the number of events

 !!    from N independent trials, each of which has

 !!    a probability of event P.

 !!      INTEGER random_binomial

 !

 !!      Method

 !

 !!     This is algorithm BTPE from:

 !

 !!         Kachitvichyanukul, V. and Schmeiser, B. W.

 !!         Binomial Random Variate Generation.

 !!         Communications of the ACM, 31, 2 (February, 1988) 216.

 !

 !!**********************************************************************

 !

 !!*****DETERMINE APPROPRIATE ALGORITHM AND WHETHER SETUP IS NECESSARY

 !

 !!     ..

 !!     .. Scalar Arguments ..

 !REAL, INTENT(IN)    :: pp

 !INTEGER, INTENT(IN) :: n

 !LOGICAL, INTENT(IN) :: first

 !INTEGER :: ival

 !!     ..

 !!     .. Local Scalars ..

 !REAL:: alv, amaxp, f, f1, f2, u, v, w, w2, x, x1, x2, ynorm, z, z2

 !REAL, PARAMETER :: ZZERO = 0.0, ZHALF = 0.5, ZONE = 1.0

 !INTEGER         :: i, ix, ix1, k, mp

 !INTEGER, SAVE   :: m

 !REAL, SAVE      :: p, q, xnp, ffm, fm, xnpq, p1, xm, xl, xr, c, al, xll,  &

 !       xlr, p2, p3, p4, qn, r, g

 !

 !!     ..

 !!     .. Executable Statements ..

 !

 !!*****SETUP, PERFORM ONLY WHEN PARAMETERS CHANGE

 !

 !IF (first) THEN

 !  p = MIN(pp, ZONE-pp)

 !  q = ZONE - p

 !  xnp = n * p

 !END IF

 !

 !IF (xnp > 30.) THEN

 !  IF (first) THEN

 !    ffm = xnp + p

 !    m = ffm

 !    fm = m

 !    xnpq = xnp * q

 !    p1 = INT(2.195*SQRT(xnpq) - 4.6*q) + ZHALF

 !    xm = fm + ZHALF

 !    xl = xm - p1

 !    xr = xm + p1

 !    c = 0.134 + 20.5 / (15.3 + fm)

 !    al = (ffm-xl) / (ffm - xl*p)

 !    xll = al * (ZONE + ZHALF*al)

 !    al = (xr - ffm) / (xr*q)

 !    xlr = al * (ZONE + ZHALF*al)

 !    p2 = p1 * (ZONE + c + c)

 !    p3 = p2 + c / xll

 !    p4 = p3 + c / xlr

 !  END IF

 !

 !!*****GENERATE VARIATE, Binomial mean at least 30.

 !

 !  20 CALL RANDOM_NUMBER(u)

 !  u = u * p4

 !  CALL RANDOM_NUMBER(v)

 !

 !!     TRIANGULAR REGION

 !

 !  IF (u <= p1) THEN

 !    ix = xm - p1 * v + u

 !    GO TO 110

 !  END IF

 !

 !!     PARALLELOGRAM REGION

 !

 !  IF (u <= p2) THEN

 !    x = xl + (u-p1) / c

 !    v = v * c + ZONE - ABS(xm-x) / p1

 !    IF (v > ZONE .OR. v <= ZZERO) GO TO 20

 !    ix = x

 !  ELSE

 !

 !!     LEFT TAIL

 !

 !    IF (u <= p3) THEN

 !      ix = xl + LOG(v) / xll

 !      IF (ix < 0) GO TO 20

 !      v = v * (u-p2) * xll

 !    ELSE

 !

 !!     RIGHT TAIL

 !

 !      ix = xr - LOG(v) / xlr

 !      IF (ix > n) GO TO 20

 !      v = v * (u-p3) * xlr

 !    END IF

 !  END IF

 !

 !!*****DETERMINE APPROPRIATE WAY TO PERFORM ACCEPT/REJECT TEST

 !

 !  k = ABS(ix-m)

 !  IF (k <= 20 .OR. k >= xnpq/2-1) THEN

 !

 !!     EXPLICIT EVALUATION

 !

 !    f = ZONE

 !    r = p / q

 !    g = (n+1) * r

 !    IF (m < ix) THEN

 !      mp = m + 1

 !      DO i = mp, ix

 !        f = f * (g/i-r)

 !      END DO

 !

 !    ELSE IF (m > ix) THEN

 !      ix1 = ix + 1

 !      DO i = ix1, m

 !        f = f / (g/i-r)

 !      END DO

 !    END IF

 !

 !    IF (v > f) THEN

 !      GO TO 20

 !    ELSE

 !      GO TO 110

 !    END IF

 !  END IF

 !

 !!     SQUEEZING USING UPPER AND LOWER BOUNDS ON LOG(F(X))

 !

 !  amaxp = (k/xnpq) * ((k*(k/3. + .625) + .1666666666666)/xnpq + ZHALF)

 !  ynorm = -k * k / (2.*xnpq)

 !  alv = LOG(v)

 !  IF (alv<ynorm - amaxp) GO TO 110

 !  IF (alv>ynorm + amaxp) GO TO 20

 !

 !!     STIRLING'S (actually de Moivre's) FORMULA TO MACHINE ACCURACY FOR

 !!     THE FINAL ACCEPTANCE/REJECTION TEST

 !

 !  x1 = ix + 1

 !  f1 = fm + ZONE

 !  z = n + 1 - fm

 !  w = n - ix + ZONE

 !  z2 = z * z

 !  x2 = x1 * x1

 !  f2 = f1 * f1

 !  w2 = w * w

 !  IF (alv - (xm*LOG(f1/x1) + (n-m+ZHALF)*LOG(z/w) + (ix-m)*LOG(w*p/(x1*q)) +    &

 !      (13860.-(462.-(132.-(99.-140./f2)/f2)/f2)/f2)/f1/166320. +   &

 !      (13860.-(462.-(132.-(99.-140./z2)/z2)/z2)/z2)/z/166320. +    &

 !      (13860.-(462.-(132.-(99.-140./x2)/x2)/x2)/x2)/x1/166320. +   &

 !      (13860.-(462.-(132.-(99.-140./w2)/w2)/w2)/w2)/w/166320.) > ZZERO) THEN

 !    GO TO 20

 !  ELSE

 !    GO TO 110

 !  END IF

 !

 !ELSE

 !!     INVERSE CDF LOGIC FOR MEAN LESS THAN 30

 !  IF (first) THEN

 !    qn = q ** n

 !    r = p / q

 !    g = r * (n+1)

 !  END IF

 !

 !  90 ix = 0

 !  f = qn

 !  CALL RANDOM_NUMBER(u)

 !  100 IF (u >= f) THEN

 !    IF (ix > 110) GO TO 90

 !    u = u - f

 !    ix = ix + 1

 !    f = f * (g/ix - r)

 !    GO TO 100

 !  END IF

 !END IF

 !

 !110 IF (pp > ZHALF) ix = n - ix

 !ival = ix

 !RETURN

 !

 !END FUNCTION random_binomial2

 !

 !

 !

 !

 !FUNCTION random_neg_binomial(sk, p) RESULT(ival)

 !

 !! Adapted from Fortran 77 code from the book:

 !!     Dagpunar, J. 'Principles of random variate generation'

 !!     Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9

 !

 !! FUNCTION GENERATES A RANDOM NEGATIVE BINOMIAL VARIATE USING UNSTORED

 !! INVERSION AND/OR THE REPRODUCTIVE PROPERTY.

 !

 !!    SK = NUMBER OF FAILURES REQUIRED (Dagpunar's words!)

 !!       = the `power' parameter of the negative binomial

 !!           (0 < REAL)

 !!    P = BERNOULLI SUCCESS PROBABILITY

 !!           (0 < REAL < 1)

 !

 !! THE PARAMETER H IS SET SO THAT UNSTORED INVERSION ONLY IS USED WHEN P <= H,

 !! OTHERWISE A COMBINATION OF UNSTORED INVERSION AND

 !! THE REPRODUCTIVE PROPERTY IS USED.

 !

 !REAL, INTENT(IN)   :: sk, p

 !INTEGER:: ival

 !

 !!     Local variables

 !! THE PARAMETER ULN = -LOG(MACHINE'S SMALLEST REAL NUMBER).

 !

 !REAL, PARAMETER    :: h = 0.7

 !REAL   :: q, x, st, uln, v, r, s, y, g

 !INTEGER:: k, i, n

 !

 !IF (sk <= XZERO .OR. p <= XZERO .OR. p >= XONE) THEN

 !  WRITE(*, *) 'IMPERMISSIBLE DISTRIBUTION PARAMETER VALUES'

 !  STOP

 !END IF

 !

 !q = XONE - p

 !x = XZERO

 !st = sk

 !IF (p > h) THEN

 !  v = XONE/LOG(p)

 !  k = st

 !  DO i = 1,k

 !    DO

 !      CALL RANDOM_NUMBER(r)

 !      IF (r > XZERO) EXIT

 !    END DO

 !    n = v*LOG(r)

 !    x = x + n

 !  END DO

 !  st = st - k

 !END IF

 !

 !s = XZERO

 !uln = -LOG(XVSMALL)

 !IF (st > -uln/LOG(q)) THEN

 !  WRITE(*, *) ' P IS TOO LARGE FOR THIS VALUE OF SK'

 !  STOP

 !END IF

 !

 !y = q**st

 !g = st

 !CALL RANDOM_NUMBER(r)

 !DO

 !  IF (y > r) EXIT

 !  r = r - y

 !  s = s + XONE

 !  y = y*p*g/s

 !  g = g + XONE

 !END DO

 !

 !ival = x + s + XHALF

 !RETURN

 !END FUNCTION random_neg_binomial

 !

 !

 !

 !FUNCTION random_von_Mises(k, first) RESULT(fn_val)

 !

 !!     Algorithm VMD from:

 !!     Dagpunar, J.S. (1990) `Sampling from the von Mises distribution via a

 !!     comparison of random numbers', J. of Appl. Statist., 17, 165-168.

 !

 !!     Fortran 90 code by Alan Miller

 !!     CSIRO Division of Mathematical & Information Sciences

 !

 !!     Arguments:

 !!     k (real)        parameter of the von Mises distribution.

 !!     first (logical) set to .TRUE. the first time that the function

 !!         is called, or the first time with a new value

 !!         for k.   When first = .TRUE., the function sets

 !!         up starting values and may be very much slower.

 !

 !REAL, INTENT(IN)     :: k

 !LOGICAL, INTENT(IN)  :: first

 !REAL     :: fn_val

 !

 !!     Local variables

 !

 !INTEGER          :: j, n

 !INTEGER, SAVE    :: nk

 !REAL, PARAMETER  :: pi = 3.14159265

 !REAL, SAVE       :: p(20), theta(0:20)

 !REAL :: sump, r, th, lambda, rlast

 !REAL (IDP)        :: dk

 !

 !IF (first) THEN! Initialization, if necessary

 !  IF (k < XZERO) THEN

 !    WRITE(*, *) '** Error: argument k for random_von_Mises = ', k

 !    RETURN

 !  END IF

 !

 !  nk = k + k + XONE

 !  IF (nk > 20) THEN

 !    WRITE(*, *) '** Error: argument k for random_von_Mises = ', k

 !    RETURN

 !  END IF

 !

 !  dk = k

 !  theta(0) = XZERO

 !  IF (k > XHALF) THEN

 !

 !!     Set up array p of probabilities.

 !

 !    sump = XZERO

 !    DO j = 1, nk

 !      IF (j < nk) THEN

 !        theta(j) = ACOS(XONE - j/k)

 !      ELSE

 !        theta(nk) = pi

 !      END IF

 !

 !!     Numerical integration of e^[k.cos(x)] from theta(j-1) to theta(j)

 !

 !      CALL integral(theta(j-1), theta(j), p(j), dk)

 !      sump = sump + p(j)

 !    END DO

 !    p(1:nk) = p(1:nk) / sump

 !  ELSE

 !    p(1) = XONE

 !    theta(1) = pi

 !  END IF ! if k > 0.5

 !END IF   ! if first

 !

 !CALL RANDOM_NUMBER(r)

 !DO j = 1, nk

 !  r = r - p(j)

 !  IF (r < XZERO) EXIT

 !END DO

 !r = -r/p(j)

 !

 !DO

 !  th = theta(j-1) + r*(theta(j) - theta(j-1))

 !  lambda = k - j + XONE - k*COS(th)

 !  n = 1

 !  rlast = lambda

 !

 !  DO

 !    CALL RANDOM_NUMBER(r)

 !    IF (r > rlast) EXIT

 !    n = n + 1

 !    rlast = r

 !  END DO

 !

 !  IF (n .NE. 2*(n/2)) EXIT         ! is n even?

 !  CALL RANDOM_NUMBER(r)

 !END DO

 !

 !fn_val = SIGN(th, (r - rlast)/(XONE - rlast) - XHALF)

 !RETURN

 !END FUNCTION random_von_Mises

 !

 !

 !

 !SUBROUTINE integral(a, b, result, dk)

 !

 !!     Gaussian integration of exp(k.cosx) from a to b.

 !

 !REAL (IDP), INTENT(IN) :: dk

 !REAL, INTENT(IN)      :: a, b

 !REAL, INTENT(OUT)     :: result

 !

 !!     Local variables

 !

 !REAL (IDP)  :: xmid, range, x1, x2,&

 !  x(3) = (/0.238619186083197_IDP, 0.661209386466265_IDP, 0.932469514203152_IDP/), &

 !  w(3) = (/0.467913934572691_IDP, 0.360761573048139_IDP, 0.171324492379170_IDP/)

 !INTEGER    :: i

 !

 !xmid = (a + b)/2._IDP

 !range = (b - a)/2._IDP

 !

 !result = 0._IDP

 !DO i = 1, 3

 !  x1 = xmid + x(i)*range

 !  x2 = xmid - x(i)*range

 !  result = result + w(i)*(EXP(dk*COS(x1)) + EXP(dk*COS(x2)))

 !END DO

 !

 !result = result * range

 !RETURN

 !END SUBROUTINE integral

 !

 !

 !

 !FUNCTION random_Cauchy() RESULT(fn_val)

 !

 !!     Generate a random deviate from the standard Cauchy distribution

 !

 !REAL     :: fn_val

 !

 !!     Local variables

 !REAL     :: v(2)

 !

 !DO

 !  CALL RANDOM_NUMBER(v)

 !  v = XTWO*(v - XHALF)

 !  IF (ABS(v(2)) < XVSMALL) CYCLE   ! Test for zero

 !  IF (v(1)**2 + v(2)**2 < XONE) EXIT

 !END DO

 !fn_val = v(1) / v(2)

 !

 !RETURN

 !END FUNCTION random_Cauchy

 !

 !

 !

 !SUBROUTINE random_order(order, n)

 !

 !!     Generate a random ordering of the integers 1 ... n.

 !

 !INTEGER, INTENT(IN)  :: n

 !INTEGER, INTENT(OUT) :: order(n)

 !

 !!     Local variables

 !

 !INTEGER :: i, j, k

 !REAL    :: wk

 !

 !DO i = 1, n

 !  order(i) = i

 !END DO

 !

 !!     Starting at the end, swap the current last indicator with one

 !!     randomly chosen from those preceeding it.

 !

 !DO i = n, 2, -1

 !  CALL RANDOM_NUMBER(wk)

 !  j = 1 + i * wk

 !  IF (j < i) THEN

 !    k = order(i)

 !    order(i) = order(j)

 !    order(j) = k

 !  END IF

 !END DO

 !

 !RETURN

 !END SUBROUTINE random_order

 !

 !

 !

 !SUBROUTINE seed_random_number(iounit)

 !

 !INTEGER, INTENT(IN)  :: iounit

 !

 !! Local variables

 !

 !INTEGER  :: k

 !INTEGER, ALLOCATABLE :: seed(:)

 !

 !CALL RANDOM_SEED(SIZE=k)

 !ALLOCATE( seed(k) )

 !

 !WRITE(*, '(a, i2, a)')' Enter ', k, ' integers for random no. seeds: '

 !READ(*, *) seed

 !WRITE(iounit, '(a, (7i10))') ' Random no. seeds: ', seed

 !CALL RANDOM_SEED(PUT=seed)

 !

 !DEALLOCATE( seed )

 !

 !RETURN

 !END SUBROUTINE seed_random_number

 !

 !

 END MODULE mode_random

mode_random::random_normal
real function random_normal()
Definition: mode_random.F90:137

mode_random
Definition: mode_random.F90:5

mode_random::init_random_seed
subroutine init_random_seed()
Definition: mode_random.F90:113