random_gamma2 Function

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 <= zero .OR. s >= one) THEN
     WRITE(*, *) 'SHAPE PARAMETER VALUE OUTSIDE PERMITTED RANGE'
     STOP
  END IF
  IF (first) THEN                        ! Initialization, if necessary
     a = one - s
     p = a/(a + s*EXP(-a))
     IF (s < vsmall) THEN
        WRITE(*, *) 'SHAPE PARAMETER VALUE TOO SMALL'
        STOP
     END IF
     c = one/s
     uf = p*(vsmall/a)**s
     vr = one - vsmall
     d = a*LOG(a)
  END IF
  DO
     CALL RANDOM_NUMBER(r)
     IF (r >= vr) THEN
        CYCLE
     ELSE IF (r > p) THEN
        x = a - LOG((one - r)/(one - p))
        w = a*LOG(x)-d
     ELSE IF (r > uf) THEN
        x = a*(r/p)**c
        w = x
     ELSE
        fn_val = zero
        RETURN
     END IF
     CALL RANDOM_NUMBER(r)
     IF (one-r <= w .AND. r > zero) THEN
        IF (r*(w + one) >= one) CYCLE
        IF (-LOG(r) <= w) CYCLE
     END IF
     EXIT
  END DO
  fn_val = x
  RETURN
END FUNCTION random_gamma2